# HG changeset patch # User Rik # Date 1460306952 25200 # Node ID f29d68e24c5a9a10b6dbfcc1dcab6412bd7e4787 # Parent cf552443c10455ede48e90c9a65ef99af6167577 ode23.m, ode45.m: More reformatting of demos to Octave coding standards. * ode23.m, ode45.m: More reformatting of demos to Octave coding standards. diff -r cf552443c104 -r f29d68e24c5a scripts/ode/ode23.m --- a/scripts/ode/ode23.m Sat Apr 09 18:58:28 2016 -0400 +++ b/scripts/ode/ode23.m Sun Apr 10 09:49:12 2016 -0700 @@ -404,6 +404,33 @@ endfunction +%!demo +%! +%! ## Demonstrate convergence order for ode23 +%! tol = 1e-5 ./ 10.^[0:8]; +%! for i = 1 : numel (tol) +%! opt = odeset ("RelTol", tol(i), "AbsTol", realmin); +%! [t, y] = ode23 (@(t, y) -y, [0, 1], 1, opt); +%! h(i) = 1 / (numel (t) - 1); +%! err(i) = 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, "-b", +%! 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", "location", "northwest"); + ## We are using the "Van der Pol" implementation for all tests that are done ## for this function. ## For further tests we also define a reference solution (computed at high @@ -590,28 +617,3 @@ %! 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"); diff -r cf552443c104 -r f29d68e24c5a scripts/ode/ode45.m --- a/scripts/ode/ode45.m Sat Apr 09 18:58:28 2016 -0400 +++ b/scripts/ode/ode45.m Sun Apr 10 09:49:12 2016 -0700 @@ -389,6 +389,33 @@ endfunction +%!demo +%! +%! ## Demonstrate convergence order for ode45 +%! tol = 1e-5 ./ 10.^[0:8]; +%! for i = 1 : numel (tol) +%! opt = odeset ("RelTol", tol(i), "AbsTol", realmin); +%! [t, y] = ode45 (@(t, y) -y, [0, 1], 1, opt); +%! h(i) = 1 / (numel (t) - 1); +%! err(i) = 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, "-b" +%! 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", "location", "northwest"); + ## We are using the "Van der Pol" implementation for all tests that are done ## for this function. ## For further tests we also define a reference solution (computed at high @@ -587,29 +614,3 @@ %! 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");