How to find second derivative as output of ode45?

9 visualizzazioni (ultimi 30 giorni)
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?

Risposta accettata

Torsten
Torsten il 22 Ago 2018
Modificato: Torsten il 22 Ago 2018
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.
  10 Commenti
I G
I G il 28 Ago 2018
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
Torsten il 29 Ago 2018
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.

Accedi per commentare.

Più risposte (0)

Prodotti

Community Treasure Hunt

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

Start Hunting!

Translated by