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

1 Ansicht (letzte 30 Tage)
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~

Akzeptierte Antwort

Jan
Jan am 20 Jun. 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 Kommentare
Yuji Zhang
Yuji Zhang am 21 Jun. 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 am 24 Jun. 2013
This is a perfect question for an internet serach engine: Simply ask your favorite engine for "Runge Kutta stepsize control"

Melden Sie sich an, um zu kommentieren.

Weitere Antworten (0)

Kategorien

Mehr zu Mathematics finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by