A problem that the Ode45 solver did not finish the computations
조회 수: 16 (최근 30일)
이전 댓글 표시
Hi, i have a problem when use the Ode45 slover. Hope someone can help me, thank you in advance!!! The problem will be described step by step as follows:
I have a boundary-value problem expressed by
I use the bvp5c solver to deal with this problem. The code is given as:
function main
solinit = bvpinit([0:0.001:L],zeros(1,6));
sol = bvp5c(@Fun,@bvpf,solinit);
plot(sol.x,sol.y(3,:))
end
function dy = Fun(x,y)
f0=230;
L=0.5;
G0=6.5282; Gw1 =1.8855e+006; Gw2 =-1.3920e-012; Gw3 =6.5986e-008; Gw4 =4.4582e-012;Gw5 =6.3222e+007;
Gz1 =1.4801e+004; Gz2 =-0.3122; Gz3 =4.5925e-014; f0=230;
dy = zeros(6,1);
dy(1) = y(2);
dy(2) = Gz1*(y(1)+y(4))-Gz2*y(6)-Gz3*y(4)*y(5);
dy(3) = y(4);
dy(4) = y(5);
dy(5) = y(6);
dy(6) = Gw1*(y(2)+y(5))+Gw2*y(5)^2+Gw3*y(4)^2+Gw3*y(4)*y(1)+Gw4*y(2)*y(5)+ ...
Gw5*y(4)*y(4)*y(5)+f0/G0*sin(pi*x/L);
end
function dy = bvpf(ya,yb)
dy = [
ya(1)
yb(1)
ya(3)
yb(3)
ya(4)
yb(4)
];
end
According to the code above, the solution is obtained readily and validated by other numerical method. The solution at x and L is obtained as
sol.y(:,1)=[0; -0.0022477; 0; 0; 0.0069115; -5.607309];
sol.y(:,end)=[0; -0.0022477; 0; 0; 0.0069115; -5.607309];
Now, i use the Ode45 solver to re-compute the problem, where the initial value is set to be sol.y(:,1):
function main
options = odeset('RelTol',1e-9,'AbsTol',1e-12);
span=[0 L]
[x,y1]=ode45(@Fun,span,sol.y(:,1),options);
end
Obviously, the solution y1(:,end) should be equal to sol.y(:,end), but ode45 solver did not even finish the computation. I don't know why this happened.
댓글 수: 3
답변 (1개)
John D'Errico
2019년 10월 21일
편집: John D'Errico
2019년 10월 21일
The answer is (probably) simple. I need point out only a couple of lines of your code that tells me it is so. (In fact, I was almost certain of it as soon as I saw your comment that ODE did not finish the computation, without even needing to look at any code at all.)
You very likely have a stiff ODE.
Usually, this seems to happen when you have different processes acting in the system you are modeling with widely varying time scales. Many chemical and biological systems seem to fall into this class of problem.
Consider these lines of code from your problem:
Gz1 =1.4801e+004; Gz2 =-0.3122; Gz3 =4.5925e-014; f0=230;
...
dy(2) = Gz1*(y(1)+y(4))-Gz2*y(6)-Gz3*y(4)*y(5);
...
G0=6.5282; Gw1 =1.8855e+006; Gw2 =-1.3920e-012; Gw3 =6.5986e-008; Gw4 =4.4582e-012;Gw5 =6.3222e+007;
...
dy(6) = Gw1*(y(2)+y(5))+Gw2*y(5)^2+Gw3*y(4)^2+Gw3*y(4)*y(1)+Gw4*y(2)*y(5)
+Gw5*y(4)*y(4)*y(5)+f0/G0*sin(pi*x/L);
So the various coefficients are either tiny, or quite large in terms of double precision. The result is that ODE45 will find it necessary to use an incredibly tiny time step to resolve those processes, but only when they become significant.
Essentially, ODE45 just gives up when the necessary time step becomes so small that no effort can be made in double precision arithmetic. (It should have told you that.)
The answer is USUALLY to switch to a stiff solver. In MATLAB, that means to use one of the solvers ode15s or ode23s. (Any solver with a name in it that ends in s.)
댓글 수: 6
참고 항목
카테고리
Help Center 및 File Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!