Solve differential equation system with ode45

Hello there,

I don't know if this is more of a mathematical problem or a programming problem...

I have the following differential equation system:

2*F+H'=0

F²+F'*H-G²-F''=0

2*F*G+H*G'-G''=0

with initial Parameters

1. F=0, G=1, H=0
2. F=0, G=0

I want to solve this with ode45 and plot the result. It should look like in the image (by Schlichting-Gersten).

How can i solve this?

Answer 1

[J. Alex Lee]

If you want to use ode45, then you have to pose as 5 first order ODEs in place of the 3 you have, by defining new variables J=G' and K=F', and then put them all in explicit form, i.e., X' = ...

Then read up on how to use ode45

[Nico Lange]

Are there other ways than ode45? Just got that as a tip, I don't know if there are any other possibilitys?

[J. Alex Lee]

yes, since you have conditions at 2 boundaries, you can look into bvp4c and bvp5c. you can also write your own FDE or FEM code. it sounds like you might want to review numerical solutions of ODEs in general.

[Nico Lange]

Okay, after some try and error i came to this code. I am still not sure how to implement the conditions at two boundaries. I've looked up bvp4c, but I dont understand how to implement the function for the boundaries.

Boundaries are:

x=0 => F=0, G=1, H=0, P=0

x=infinite => F=0, G=0

I added a fourth differential equation to the three above P'+H*H'-H''=0.

After that I've build 6 first order ODEs, because of 6 conditions, where I eliminated F (not sure if this was right)

x(1) ' = H' = y(2)

x(2) ' = H'' = y(3)

x(3) ' = H''' = ( (-0.5*y(2))^2 - 0.5*y(3)*y(1) - y(4)^2 ) / (-0.5)

x(4) ' = G' = y(5)

x(5) ' = G'' = 2*(-0.5*y(2))*y(4) + y(1)*y(5)

x(6) ' = H'' - HH' = y(3) - y(1)*y(2)

How do I implement the conditions for the boundaries?d (I've tried it with y0 like in the code).

tol=1e-3;   
x=0;              
dx=0.1;         
Xf=4;            
tspan=(0:dx:Xf);   
Nt=Xf/dx+1;

for i=1:10000       
  x=x+0.0001;
  y0=[x;1;0;0;0;0];     
  [t,y]=ode45(@(t,y) odefcn(tspan,y),tspan,y0);
  y2=(y(Nt,2));
  if abs(y(Nt,2)-1.0)<tol, break, end
end

figure(1)
plot(t,y(:,1),t,y(:,2),t,y(:,3),t,y(:,4),t,y(:,5),t,y(:,6),'LineWidth',2)
axis([0 4 0 1]);
grid on;


function dydt = odefcn(tspan,y)
%        H    = y(1)
%x(1)' = H'   = y(2)
%x(2)' = H''  = y(3)
%x(3)' = H''' = ((-0.5H')²+(-0.5H'')H-G²)/(-0.5) = ...
%               ((-0.5*y(2))²-0.5*y(3)*y(1)-y(4))/(-0.5)
%        G    = y(4)
%x(4)' = G'   = y(5)
%x(5)' = G''  = 2*(-0.5H')G+HG' = 2*(-0.5*y(2))*y(4)+y(1)*y(5)
%x(6)' = H'' - HH' = y(3)-y(1)*y(2)
dydt = [y(2);
        y(3); 
        ((-0.5*y(2))^2-0.5*y(3)*y(1)-y(4)^2)/(-0.5);
        y(5);
        2*(-0.5*y(2))*y(4)+y(1)*y(5);
        y(3)-y(1)*y(2)];
end

[J. Alex Lee]

What is P? Are you now solving a different set of equations?

If your boundary is at infinity, one strategy using ode45 (or any other odeXX) is to guess the unknown BCs at the first boundary and integrate until the actual BCs are met (using an event function) and then check the derivatives, because they also ought to be zero.

I'm not sure how to deal with infinite boundaries with bvp4c/5c because i don't use them much...it is a common type of problem though, don't hte documentation have an example about it?

[J. Alex Lee]

You should move your "Answer" back into these comments...

Well, ok, so P is only coupled to the system one-way, so you can compute it later, so strictly it is not necessary to include. But if you want to compute it even if it's not in the example solution curves, it's probably easier just to include in the system so that you don't have to post-process later.

I agree that the link you posted is relevant to find approximate solutions that are consistent with the infinite boundary situation that you have.

The "residual form" of any condition is to rearrange the equation in the form 0 = ...That way you can pose as a root-finding operation. Using bvp4/5c you won't need to worry about the details of how to solve that. On the other hand if you want to use shooting with ode45 or some variant, you will need to solve at a least a 2-dimensional root finding problem (for the 2 missing initial conditions at x=0). Bottom line, residual form is just equation rearranged so that the LHS is 0, and you replace zero with "res", which is the thing you are asking bvp4c to set to 0.

Based on your comments and questions, I would recommend learning how to use bvp4c rather than the other approaches I mentioned.

Answer 2

[Nico Lange]

No, the full system of ODEs ist

2*F+H'=0

F²+F'*H-G²-F''=0

2*F*G+H*G'-G''=0

P'+H*H'-H''=0

with initial Parameters

x=0 => F=0, G=1, H=0, P=0

x=infinite => F=0, G=0

but I wasn't sure if i needed the last one. I just added it yesterday, just in case that it is necessary.

https://de.mathworks.com/help/matlab/math/verify-bvp-consistency-using-continuation.html

Yes there is an example, but I am not sure how to work with res = fsbc(), beacuse I never needed it and I don't know how to use the function with bvp4c. I also never worked with bvp4c before, so I am not sure how to use it correctly :'D

[James Tursa]

This is confusing to me, and as posted seems unsolvable. You have these dependent variables: F, G, H, P, and the independent variable x. The highest order derivatives appearing in your equations are F'', G'', H'', and P'. This would mean you have a 2+2+2+1 = 7th order system of equations. The state for this system would be (F,F',G,G',H,H',P). Call this Y. The differential equations governing this state would be:

dY(1)/dx = dF/dx = F' = Y(2)

dY(2)/dx = dF'/dx = F'' = (the expression from above solved for F'')

dY(3)/dx = dG/dx = G' = Y(4)

dY(4)/dx = dG'/dx = G'' = (the expression from above solved for G'')

dY(5)/dx = dH/dx = H' = Y(6)

dY(6)/dx = dH'/dx = H'' = (the expression from above solved for H'')

dY(7)/dx = dP/dx = P' = ???

You have posted nothing for that last one, P', in terms of just the state, so you don't have enough differential equations to solve this system.

Also, you list values for some of the states at x=0 and others at x=infinity. So this is a boundary value problem, not an initial value problem. ode45( ) can't be used for this since it is not an initial value problem.

[Nico Lange]

I am not sure if I need the equation with P to plot the curves shown above. I just want to plot F, G and H. I don't have more differetial equations, just the one I showed you. What if I just leave P'+H*H'-H''=0 out?

So even if I try this with bvp4c, I am not sure if I made the correct system of ODEs.

[J. Alex Lee]

oh i see now, there's an H'' in the last equation. What's missing right now is not another differential equation, but a boundary condition. the most straightforward missing BC is then a second condition on H; there's already 1 condition for P.

but if Nico doesn't care about solving for P, then I agree the last equation can be left out, and the original post should have enough conditions to solve.

Yes it is a BVP, but it can be solved using ODE45 + shooting method posing a root-finding problem in the 2 unknown BC's at x=0. But, it is perhaps less straightforward to do this than to use bvp4c/5c with the domain-stretching continuation method linked, especially if there's not much familiarity with numerical methods for ODE solving and coding.

[Nico Lange]

Okay, so I've tried it this way: I made 8 ODEs with F, H,H', H'', H''', G, G', G'' and looked up bvp4c. Code below. Problem is, that i get the Error

Error using bvp4c (line 251)
Unable to solve the collocation equations -- a singular Jacobian
encountered.

Error in Testneu (line 4)
  sol = bvp4c(@fsode, @fsbc, solinit);

What do I do wrong? The ODEs should be correct, I've tried it with ode45 before and I've got plots similar to the picture above, but still wrong.

  clc; close all; clear all;

infinity = 4;
solinit = bvpinit(linspace(0,infinity,5),[0 -1.2 2 -0.5 0 0 0 1]);

 sol = bvp4c(@fsode, @fsbc, solinit);
 x = sol.x;
 f = sol.y;

figure(1)
%plot(x,f(:,6),x,f(:,8),x,f(:,7),'LineWidth',2)
axis([0 4 0 1]);
grid on;
legend('F','G','-H');
xlabel('\zeta=z\surd(\omega/\nu)');
ylabel('H,G,F');

function dfdeta = fsode(x,f)
dfdeta = [f(2);
        f(3);
        ((-0.5*f(2))^2-0.5*f(3)*f(1)-f(4)^2)/(-0.5);
        f(5);
        2*(-0.5*f(2))*f(4)+f(1)*f(5);
        -0.5*f(2);
        -f(1);
        f(4)];      
end

function res = fsbc(f0,finf)
res =[f0(1)-0;
      finf(1)-0.88446;
      f0(2)-0;
      finf(2)-0;
      f0(3)-0;
      f0(4)-1;
      finf(4)-0;
      f0(5)-0];
end

[J. Alex Lee]

it does you no good to add more odes... stick with the 5 ode's back when you didn't need the P and follow the examples of bvp4c/5c.