How to find second derivative as output of ode45?

I am using ode45 function to find numerical solution for my system of equations, where I have 4 equations and 4 variables, with command:

sol=ode45(@fun,[1 0],[1; 0; 0; 0])

where time span is going from 1 to 0, and initial conditions are 1,0,0,0.

So, I need to find numerically the first and the second derivative of my values according to time. For the first derivative I have verified code

[~,SXINT]=deval(sol,sol.x);

where in SXINT are stored first derivatives of all 4 variables.

But I don`t know how to find the second derivative numerically. I tried:

[~,SXINT2]=deval(SXINT,sol.x);

but I got error:

SXINT must be a structure returned by a differential equation solver.

How to make that SXINT be in that structure, or is there some other way to find the second derivative, but with the same precision as deval method is with ode45 solver?

Answer 1

[Torsten]

To get the second derivative, the first derivative must be part of the solution "sol".

So you have to solve 8=4*2 equations, namely your original system and new algebraic equations given by

y(5)-dy(1) = 0
y(6)-dy(2) = 0
y(7)-dy(3) = 0
y(8)-dy(4) = 0

This way, "sol" has stored the first-order derivatives of y(1)-y(4) in y(5)-y(8) which can now be differentiated by a call to "deval".

Best wishes

Torsten.

[I G]

It works. But now in sol I have solution of equations and its first derivatives and I don't know how to extract them. I need only the first four solution from sol to get graphics of solution of ode45, or only the second four solutions (it means first derivatives) and than to use it to obtain the second derivatives?

Exactly what to call on ???

[~,SXINT]=deval(???,sol.x);

I have tried with sol(1), it doesn't work.

[I G]

Ok, I got it, but I have the problem because as a result I got that my terms in matrix are the same, that means the first and the fifth columns are the same. I don't know how, here is the code:

sol=ode45(@fun,[1 0],[1; 0; 0; 0; 0; 0; 0; 0;])
[~,SXINT]=deval(sol,sol.x);
dp=SXINT(1,:);
dp2=SXINT(5,:)

plot(sol.x,do,sol.x,dp1,'x');

According to code, as result of deval I should have the first derivatives in the first four columns, and the fourth derivatives in the columns 5-8. Why are they the same?

In my fun.m, according to you suggestions I added:

f(5)=f(1);
f(6)=f(2);
f(7)=f(3);
f(8)=f(4);

*I have the same result when I use only the second four terms to get derivatives:

[~,SXINT]=deval(sol,sol.x,[5 6 7 8]);

Maybe it solves and terms f(5)...f(8) in the same way, as the first four terms. But how to tell it, that I need that values of derivatives when I get out of function, and not values of variables?

[Torsten]

Change the call to ODE45 as:

M = diag([1 1 1 1 0 0 0 0]);
options = odeset('Mass',M);
sol = ode45(@fun,[1 0],[1; 0; 0; 0; 0; 0; 0; 0;],options)

Change fun.m (last four lines) as

f(5)=y(5)-f(1);
f(6)=y(6)-f(2);
f(7)=y(7)-f(3);
f(8)=y(8)-f(4);

Then SXINT(1,:) = y1', ...,SXINT(4,:)=y4', SXINT(5,:)=y1'',...,SXINT(8,:)=y4''.

Best wishes

Torsten.

[I G]

I got NaN values for the all second derivatives.

[I G]

Where should I look for the problem?

[Torsten]

Use

M = diag([1 1 1 1 0 0 0 0]);
options = odeset('Mass',M,'MassSingular','yes');
sol = ode15s(@fun,[1 0],[1; 0; 0; 0; 0; 0; 0; 0],options)

Best wishes

Torsten.

[I G]

I am using ode45, not ode15s, because I use Runge Kutta method. From that point of view will ode15s be of the same precision as ode 45?

This is results of you suggestions:

- with ode15s, if I compare results in SXINT(4,:) and SINT(8,:) they are the same, so that is not good, because I am looking for the first and the second derivatives-they not need to be the same

- with ode15s, if I compare results in SXINT(4,:) and SXINT(8,:) they are not the same, but I am not sure is that thing which I am looking for, and it is ode15s, I need ode45

- if I use ode45, it is

MassSingular cannot be 'yes' for this solver.

This is my function file, where p is variable, f is the first derivative of p:

function f=fun_z1(z,p)

beta=1;
f=zeros(8,1);
ri=2;

R=ri-z*(ri-1);
f(1)=-32*beta/(R.^4*p(1));
f(2)=(-8*f(1)/R-f(1)*p(2))/p(1);
f(3)=(-p(2)*f(2)+(-8*f(2)/R-8*f(1)/(R.*R*p(1)))-f(1)*p(3))/p(1);
f(4)=(-f(2)*p(3)-f(3)*p(2)+8*(-f(3)/R-f(2)/p(1)+p(2)*f(1)/(p(1).*p(1))) -f(1)*p(4))/p(1);
f(5)=p(5)-f(1);
f(6)=p(6)-f(2);
f(7)=p(7)-f(3);
f(8)=p(8)-f(4);

This is code for comparing the first derivative and potential second derivative:

M=diag([1 1 1 1 0 0 0 0]);
options=odeset('Mass',M,'MassSingular','yes');

sol=ode15s(@fun_z1,[1 0],[1; 0; 0; 0;0;0;0;0],options);

[SINT,SXINT] = deval(sol,sol.x);

dp3dz_num_sol=SXINT(4,:);
dp32dz_num_sol=SINT(8,:);

plot(sol.x,dp3dz_num_sol,sol.x,dp32dz_num_sol,'x')

What do you think, does ode15s is ok to use instead ode45, if I am looking for Runge Kutta solver?

[Torsten]

Use ODE15S ; in most cases, it performs better than ODE45.

For completeness, I'll list the contents of SINT and SXINT from the call

[SINT,SXINT] = deval(sol,sol.x);

:

SINT(1,:) : p(1)
SINT(2,:) : p(2)
SINT(3,:) : p(3)
SINT(4,:) : p(4)
SINT(5,:) : d/dz p(1)
SINT(6,:) : d/dz p(2)
SINT(7,:) : d/dz p(3)
SINT(8,:) : d/dz p(4)
SXINT(1,:) : d/dz p(1)
SXINT(2,:) : d/dz p(2)
SXINT(3,:) : d/dz p(3)
SXINT(4,:) : d/dz p(4)
SXINT(5,:) : d^2/dz^2 p(1)
SXINT(6,:) : d^2/dz^2 p(2)
SXINT(7,:) : d^2/dz^2 p(3)
SXINT(8,:) : d^2/dz^2 p(4)

Best wishes

Torsten.

[I G]

For this case I also have analytical solution, and my goal is to compare analytical solution and numerical solution of the second derivatives. But I got in the first case - first equation in my system -that analytical and numerical solution are not the same. It means they are similar, but there exist some difference which is not negligible. Here is my code (attached "za_sajt.m"), that would be "1. aprox" plot.

I have tried with odeset options 'RelTol' and 'AbsTol' but I got just worse matching, even not the same trend with this options. Do you have some advice?

[Torsten]

The way I usually calculate higher derivatives of ODE variables is simply differentiating the right-hand side of the ODE equations with respect to the independent variable:

If

dyi/dz = fi(z,y1(z),y2(z),...,yn(z)),

then

d^2yi/dz^2 = partial(fi)/partial(z) + sum_j (partial(fi)/partial(yj)*dyj/dz)

E.g. for your first equation

beta = 1.0;
ri = 0.3;
R = ri-z*(ri-1);
dR/dz = -ri;
dp1/dz = -32*beta/(R^4*p(1));
d^2p1/dz^2 = 32*beta*(4*R^3*dRdz*p(1)+R^4*dp1/dz)/(R^4*p(1))^2

Best wishes

Torsten.