A problem that the Ode45 solver did not finish the computations

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.

[David Wilson]

Did you forget the three dots (line continuation) in the ODE function ?

 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);

In any case, the solving the ODE could be (numerically) unstable, whereas solving the BVP is actually OK.

[Xuan Ling Zhang]

Dear Wilson,

sorry, I just forget the three dots in the above edit, and i'm sure that there is no problem on the ODE function in my Matlab m file.

[Xuan Ling Zhang]

Dear Wilson,

Thanks for your good suggestion in the comments below, according to which i will try it.

Best wishes.

Yours sincerely,

Answer 1

[John D'Errico]

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.)

https://en.wikipedia.org/wiki/Stiff_equation

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.)

[Xuan Ling Zhang]

Dear D'Errico,

i had tried the ode15s solver before, now i have just tried the ode23s solver according to your suggestion. but it still doesn't work.

[Xuan Ling Zhang]

Actually, this equations are used to describe a beam's nonlinear problem，which have tiny coefficients. Now i give up the nonlinear terms to investigate the linear problem. that is keep the coefficients Gz1,Gz2 and Gw1, which can be seen that the equations are not stiff equations. but the similar result emerge as follows：

[David Wilson]

I'm not so sure stiffness is necessarily the main issue here (although it could be a complicating factor). I too tried various solvers, and the key point is not that the ode45 did not finish, but that it diverged & threw an exception, which is quite a different problem.

I did not fully investigate this, but I suspect that the problem is similar to a linear problem where the homogeneous solution has positive (unstable) eigenvalues, but the particular integral nullifies them making the actual solution technically stable, but any error will be unstable. Consequently any marching scheme (such as odexxx) will fail. Collocation could well work.

I suspect that you could integrate the system backwards in time from the known final conditions which may then work. (Use a variable substitution tau = -t, and re-write your ODES in tau.) Or you could try using the bvp4c in "ODE" mode, i.e. specify your 6 boundary conditions all at the same boundary and re-compute.

[Steven Lord]

When I try setting one additional option to see what the solution looks like:

options = odeset('RelTol',1e-9,'AbsTol',1e-12, 'OutputFcn', @odeplot);

I see your function falls off a cliff at around 0.02655434. By the time I stopped it at x = 0.0265543425513827, the value of the third component of the solution was around -0.00688 and the value of the sixth was -3.97e34.

[Xuan Ling Zhang]

Dear Lord,

Yes, this is what happened that i don't know why. Do you have some suggestions？

[Xuan Ling Zhang]

Dear Wilson, the first scheme of 'integrate the system backwards' that i have tried but fails.