Mercurial > octave
changeset 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 | 3933d5073baa |
| children | fe1447ae68cf |
| files | scripts/ode/ode23.m scripts/ode/ode45.m |
| diffstat | 2 files changed, 51 insertions(+), 0 deletions(-) [+] |
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");
--- a/scripts/ode/ode45.m Fri Apr 08 16:51:05 2016 -0400 +++ b/scripts/ode/ode45.m Sat Apr 09 06:36:53 2016 +0200 @@ -587,3 +587,29 @@ %! opt = struct ("foo", 1); %! [t, y] = ode45 (@fpol, [0 2], [2 0], opt); + +%!demo +%! +%! # Demonstrate convergence order for ode45 +%! tol = 1e-5 ./ 10.^[0:8]; +%! for ii = 1 : numel (tol) +%! opt = odeset ("RelTol", tol(ii), "AbsTol", realmin); +%! [t, y] = ode45 (@(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)) .^ 4 .* tol(end), "k--", +%! h, (h/h(end)) .^ 5 .* tol(end), "k-") +%! axis tight +%! xlabel ("h"); +%! ylabel ("err(h)"); +%! title ("Convergence plot for ode45"); +%! legend ("imposed tolerance", "ode45 (relative) error", +%! "order 4", "order 5");