Solve differential equation system with ode45

조회 수: 2 (최근 30일)
Nico Lange
Nico Lange 2020년 12월 28일
댓글: J. Alex Lee 2021년 1월 3일
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?

답변 (2개)

J. Alex Lee
J. Alex Lee 2020년 12월 28일
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
J. Alex Lee
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
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.

댓글을 달려면 로그인하십시오.


Nico Lange
Nico Lange 2020년 12월 30일
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
Nico Lange
Nico Lange 2021년 1월 2일
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
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.

댓글을 달려면 로그인하십시오.

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by