manually run ode45 step by step - and impose filter to the half-way solution

1 vue (au cours des 30 derniers jours)
Yuji Zhang
Yuji Zhang le 20 Juin 2013
Hi everybody~
I'm solving a ode:
%%code 1
t = [tt1 tt2 tt3 ... tt100];
[t, y] = ode45(@equation, t, y1, options); % y(t=1)=y1
Now I wanna to run it step by step. I wanna verify the above code is equivalent to
%%code 2
t = [tt1 tt2 tt3 ... tt100];
for i = 1:99
[t(i+1) y(i+1)] = ode45(@equation, [t(i) t(i+1)], y(i), options);
end
The reason I want step-by-step is, I want to impose filters to y in the middle. Say, at t=t(5), I do y(5) = y(5)*filter, and then continue to solve y(6).
Let me know whether code 1 is equivalent to code 2.
Any thoughts of other ways to do this work is appreciated. I think this is called "numerical filtering" - but didn't find anything to read about it. Any recommended reading is appreciated. if true % code endThanks~

Réponse acceptée

Jan
Jan le 20 Juin 2013
The calculations are equivalent, but not exactly identical: While in the first version ODE45 uses the stepsize from the former step after reaching one of the intermediate points, the 2nd method restarts the integrations with a new estimated stepsize. Due to rounding and local discretization errors, the results wil be slightly different - except if the solution is instable and the small deviations are amplified, which can cause large differences. But then the "result" of the integration is doubtful at all.
  2 commentaires
Yuji Zhang
Yuji Zhang le 21 Juin 2013
Hi Jan~
Thank u for your help.
I learn from "help" that, t= t(i) is just the "record points" I want y(i) to be recorded. code45 automatically determines internal step-length of delta_t. You were talking about how the step-length is determined. Is there something I can read about it?
Let me know. Thank u!
Jan
Jan le 24 Juin 2013
This is a perfect question for an internet serach engine: Simply ask your favorite engine for "Runge Kutta stepsize control"

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Mathematics dans Help Center et File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by