Solve differential equation system with ode45

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

Nico Lange 2020-12-29

在 MATLAB Online 中打开

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 2020-12-29

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 2020-12-30

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 2020-12-30

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/703822-solve-differential-equation-system-with-ode45#answer_587578

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

5 个评论
显示 3更早的评论隐藏 3更早的评论

James Tursa 2020-12-30

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 2021-1-2

在 MATLAB Online 中打开

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 2021-1-3

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.

请先登录，再进行评论。

Solve differential equation system with ode45

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

5 个评论
显示 3更早的评论隐藏 3更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

Community Treasure Hunt

Solve differential equation system with ode45

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（2 个）

5 个评论 显示 3更早的评论隐藏 3更早的评论

5 个评论 显示 3更早的评论隐藏 3更早的评论

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论

5 个评论
显示 3更早的评论隐藏 3更早的评论