MATLAB Answers

Time derivative of parameters within ODE solvers

3 views (last 30 days)
I have an ODE which has a parameter whose 1st and 2nd order time derivatives are also included in the ODE:
function dy = ODE(t,y)
f1 = myfun(y(t));
dy = y + f1 + df1/dt + ddf1/dt^2;
end
Unfortunately, the function of analytical derivatives of myfun is not available. Therefore, df1 and ddf1 can be computed numerically, only. Given that the time step in Matlab ode solvers is not fixed, I wonder if there is a way to numerically compute df1 and ddf1.

  0 Comments

Sign in to comment.

Accepted Answer

Torsten
Torsten on 23 Nov 2018
Edited: Torsten on 23 Nov 2018
function dy = ODE(t,y)
dt = 1e-8;
fm = myfun(t-dt);
f = myfun(t);
fp = myfun(t+dt);
df = (fp - fm) / (2 * dt);
ddf = (fp - 2 * f + fm) / dt^2;
dy = y + f + df + ddf;
end

  6 Comments

Show 3 older comments
Torsten
Torsten on 23 Nov 2018
Does "myfun" depend explicitly on t or only via y ? I.e. can myfun look like y^2+t^2 or only like y^2 ?
Torsten
Torsten on 23 Nov 2018
Let
z = f(y)
the value that "myfun" returns for argument y.
Then
dz/dt = df/dy * dy/dt
d^2z/dt^2 = d^2f/dy^2 * (dy/dt)^2 + df/dy * d^2y/dt^2
Inserting into your differential equation gives
dy/dt = y + f + df/dy * dy/dt + d^2f/dy^2 * (dy/dt)^2 + df/dy * d^2y/dt^2
or
df/dy * d^2y/dt^2 + (df/dy - 1) * dy/dt + d^2f/dy^2 * (dy/dt)^2 + y + f = 0
Now you can approximate df/dy and d^2f/dy^2 as I suggested above and solve the system (z1 = y, z2 = dy/dt)
z1' = z2
z2' = -((df/dy - 1) * z2 + d^2f/dy^2 * (z2)^2 + z1 + f)/(df/dy)
using ODE45, e.g.

Sign in to comment.

More Answers (0)

Tags

Products


Release

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by