Calculates the solution of an ODE $ \dot x(t) = f(t,x) $ for a given Butcher parametrisation.
| func | The function $f$ defining the ODE $ \dot x(t) = f(t,x). $ |
| t0 | The start of the evaluation interval |
| t1 | The end of the evaluation interval |
| dt | Intitial iteration step (might be adjusted later, c.f. above). |
| butcher | The butcher tables to parametrise the Runge-Kutta method. |
| adapt_set | Setting of the step adaptive steps. |
| unit_hint | The method tries to guess all the right quantities from the parameters, but zero might trick it. You can give here the right units [t_units, [x0_unit, x1_unit, ...]] |
| return value | Either the value of the solution in format [t0, x00, x10, ...; t1, x01, x11, ...], or the coefficents for the interpolation when interpolation is calculated. |
Our example in this documentation will be the autmonomous $ \ddot x=x $.
Fix step iteration.
In this case, the iteration goes through with the predefined t-step between $t_0$ and $t_1$. The stucture 'butcher' requires only a 'matrix', a 'steps', and a 'coeff' component. They are coefficients of the previous function evaluations for the new x in the new evaluation, the time steps for the new evaluations, and the final coefficients for the eventual step respectively. See the return values of 'odeRungeKutta3()' and 'odeRungeKutta4()' as examples.
soln = odeRK(@(t,x) [x[1],-x[0]/1s^2], 0, 5s, .01s, [1e-21m, 1m/s], odeRungeKutta4());
soln.size()
plot(soln[...;0]/1s, soln[...;1]/1m)
Often times one needs an error control on the t-steps. There are embedded RK methods which estimate errors based on a different order method derived from the same function evaluation matrix. The parameter 'butcher' will require the corresponding final coefficents in the component 'corr', and also a 'corr_fn' which returns the error estimation from the last two (hence allowing simple control system methods) difference vector of the two different order results. See the return values of 'odeDormandPrince54()' and 'odeDormandPrince853' as examples.
To enable the adaptive version, one need to define the parameter 'adapt_set'. It takes the following components:
soln = odeRK(@(t,x) [x[1],-x[0]/1s^2], 0, 5s, .01s, [1e-21m, 1m/s], odeDormandPrince853(), {});
soln.size()
plot(soln[...;0]/1s, soln[...;1]/1m)
The clear advantage of high-order adaptive methods are that they can make huge steps, wherever it is possible. Sometimes one needs the values inbeetween the iteration steps. Standard interpolation techniques apply here too. The stucture 'butcher' will need additional fields 'interp', 'ip_coeff', 'ip_steps'. They will correspond to the 'matrix', 'coeff', and 'steps' of the base method. The first four coefficients of the interpolation will be calculated from the iteration step, and the additional coefficients will be the evaluations of 'ip_coeff' for the Horner like polynomial evaluation: $$ p(t) = c_0 + t(c_1 + (1-t)*(c_2 + t*(c_3 + \dots))), $$ where $ t $ is the fraction of the interpolation step. This evaluation is performed by the function 'odeRKEvalInterp()'.
The method is switched on by setting the component 'interp' of the parameter 'adapt_set' true.
f = @(t,x) [x[1],-x[0]/1s^2];
t0 = 0;
t1 = 5s;
dt = .01s;
x0 = [1e-21m, 1m/s];
B = almafa.odeDormandPrince853();
sparse = almafa.odeRK(f, t0, t1, dt, x0, B, {});
interp = almafa.odeRK(f, t0, t1, dt, x0, B, {interp:true});
j0 = 21; j1 = 29;
t = sparse[j0;0]..1ms..sparse[j1;0]
x = (@(_) almafa.odeRKEvalInterp(interp,_)[0])(#t);
plot(sparse[j0..j1;0]/1s,sparse[j0..j1;1]/1m).curve(t/1s,x/1m)