Mercurial > octave
changeset 20847:ddc18b909ec7
codesprint: ode: first and second order interpolator for dense_output
*scripts/ode/private/integrate_adaptive: integrator passes
the function to the interpolator in case it needs
further function evaluations
*scripts/ode/private/runge_kutta_interpolate:
use 'interp1 (.., 'linear') for order 1 and
hermite_quadratic_interpolotion for order 2
*scripts/ode/private/starting_stepsize: correctly handle
the case in wich func returns a cell
| author | jcorno <jacopo.corno@gmail.com> |
|---|---|
| date | Fri, 20 Nov 2015 15:38:10 +0100 |
| parents | 0bb9b4e99ea0 |
| children | 56d36905893f |
| files | scripts/ode/private/integrate_adaptive.m scripts/ode/private/runge_kutta_interpolate.m scripts/ode/private/starting_stepsize.m |
| diffstat | 3 files changed, 29 insertions(+), 8 deletions(-) [+] |
line wrap: on
line diff
--- a/scripts/ode/private/integrate_adaptive.m Sat Dec 12 15:19:13 2015 +0100 +++ b/scripts/ode/private/integrate_adaptive.m Fri Nov 20 15:38:10 2015 +0100 @@ -150,8 +150,8 @@ iout = max (t_caught); x(:, t_caught) = ... runge_kutta_interpolate (order, [t_old t_new], [x_old x_new], - tspan(t_caught), new_k_vals, dt, - options.funarguments{:}); + tspan(t_caught), new_k_vals, dt, func, + options.funarguments); istep += 1;
--- a/scripts/ode/private/runge_kutta_interpolate.m Sat Dec 12 15:19:13 2015 +0100 +++ b/scripts/ode/private/runge_kutta_interpolate.m Fri Nov 20 15:38:10 2015 +0100 @@ -17,14 +17,15 @@ ## along with Octave; see the file COPYING. If not, see ## <http://www.gnu.org/licenses/>. -function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals, dt, args) +function u_interp = runge_kutta_interpolate (order, z, u, t, k_vals, dt, func, args) switch (order) ## FIXME: Can interpolations for orders 1-4 simply be deleted? 2015-10-14. - #{ + case 1 - u_interp = linear_interpolation (z, u, t); + u_interp = interp1 (z, u', t, 'linear'); + case 2 if (! isempty (k_vals)) der = k_vals(:,1); @@ -32,6 +33,8 @@ der = feval (func, z(1) , u(:,1), args); endif u_interp = quadratic_interpolation (z, u, der, t); + + #{ case 3 u_interp = ... hermite_cubic_interpolation (z, u, k_vals, t); @@ -42,7 +45,6 @@ [u(:,i-1) u(:,i)], k_vals, tspan(counter)); #} - case 5 ## ode45 with Dormand-Prince scheme: u_interp = hermite_quartic_interpolation (z, u, k_vals, t); @@ -69,6 +71,19 @@ endfunction +## The function below can be used in an ODE solver to interpolate the solution +## at the time t_out using 2th order hermite interpolation. +function x_out = quadratic_interpolation (t, x, der, t_out) + + # coefficients of the parabola + a = -(x(:,1) - x(:,2) - der(:).*(t(1)-t(2))) / (t(1) - t(2))^2; + b = der(:) - 2*t(1).*a; + c = x(:,1) - a*t(1)^2 - b*t(1); + + # evauate in t_out + x_out = a*t_out.^2 + b*t_out + c; + +endfunction ## The function below can be used in an ODE solver to interpolate the solution ## at the time t_out using 4th order hermite interpolation.
--- a/scripts/ode/private/starting_stepsize.m Sat Dec 12 15:19:13 2015 +0100 +++ b/scripts/ode/private/starting_stepsize.m Fri Nov 20 15:38:10 2015 +0100 @@ -43,6 +43,9 @@ ## compute norm of the function evaluated at initial conditions y = func (t0, x0); + if (iscell (y)) + y = y{1}; + endif d1 = AbsRel_Norm (y, y, AbsTol, RelTol, normcontrol); if (d0 < 1e-5 || d1 < 1e-5) @@ -55,9 +58,12 @@ x1 = x0 + h0 * y; ## approximate the derivative norm + yh = func (t0+h0, x1); + if (iscell (yh)) + yh = yh{1}; + endif d2 = (1 / h0) * ... - AbsRel_Norm (func (t0+h0, x1) - y, - func (t0+h0, x1) - y, AbsTol, RelTol, normcontrol); + AbsRel_Norm (yh - y, yh - y, AbsTol, RelTol, normcontrol); if (max (d1, d2) <= 1e-15) h1 = max (1e-6, h0*1e-3);