Time derivative of parameters within ODE solvers
4 views (last 30 days)
Show older comments
Seyed Ali Baradaran Birjandi
on 22 Nov 2018
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
Accepted Answer
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
Torsten
on 23 Nov 2018
Edited: 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.
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!