Mercurial > octave
diff scripts/ode/ode23.m @ 21596:07d30e6fcfde
Add demos to ode23 and ode45
* scripts/ode/ode{23,45}.m : Add new demo demonstrating
order of convergence.
author | Carlo de Falco <carlo.defalco@polimi.it> |
---|---|
date | Sat, 09 Apr 2016 06:36:53 +0200 |
parents | ecce63c99c3f |
children | f29d68e24c5a |
line wrap: on
line diff
--- a/scripts/ode/ode23.m Fri Apr 08 16:51:05 2016 -0400 +++ b/scripts/ode/ode23.m Sat Apr 09 06:36:53 2016 +0200 @@ -590,3 +590,28 @@ %! opt = struct ("foo", 1); %! [t, y] = ode23 (@fpol, [0 2], [2 0], opt); +%!demo +%! +%! # Demonstrate convergence order for ode23 +%! tol = 1e-5 ./ 10.^[0:8]; +%! for ii = 1 : numel (tol) +%! opt = odeset ("RelTol", tol(ii), "AbsTol", realmin); +%! [t, y] = ode23 (@(t, y) -y, [0, 1], 1, opt); +%! h(ii) = 1 / (numel (t) - 1); +%! err(ii) = norm (y .* exp (t) - 1, inf); +%! endfor +%! +%! # Estimate order numerically +%! p = diff (log (err)) ./ diff (log (h)) +%! +%! %! # Estimate order visually +%! figure; +%! loglog (h, tol, '-ob', h, err, +%! h, (h/h(end)) .^ 2 .* tol(end), "k--", +%! h, (h/h(end)) .^ 3 .* tol(end), "k-") +%! axis tight +%! xlabel ("h"); +%! ylabel ("err(h)"); +%! title ("Convergence plot for ode23"); +%! legend ("imposed tolerance", "ode23 (relative) error", +%! "order 2", "order 3");