Mercurial > octave
diff scripts/ode/ode15s.m @ 29456:ebc3f80673f0
ode15s.m: Fix error if Jacobian in options is a constant matrix (bug #60240).
* scripts/ode/ode15s.m: If the mass matrix is empty, put the identity for
df/dy' into options.Jacobian. Add a linear system to the BISTs.
| author | Markus Meisinger <chloros2@gmx.de> |
|---|---|
| date | Sun, 21 Mar 2021 10:53:11 +0100 |
| parents | 7854d5752dd2 |
| children | 363fb10055df |
line wrap: on
line diff
--- a/scripts/ode/ode15s.m Sat Mar 20 15:43:04 2021 +0100 +++ b/scripts/ode/ode15s.m Sun Mar 21 10:53:11 2021 +0100 @@ -257,7 +257,7 @@ endif ## Use sparse methods only if all matrices are sparse - if (! options.havemasssparse) + if (! isempty (options.Mass)) && (! options.havemasssparse) options.havejacsparse = false; endif @@ -271,7 +271,15 @@ options.havejacfun); options.havejacfun = true; else # All matrices are constant - options.Jacobian = {[- options.Jacobian], [options.Mass]}; + if (! isempty (options.Mass)) + options.Jacobian = {[- options.Jacobian], [options.Mass]}; + else + if (options.havejacsparse) + options.Jacobian = {[- options.Jacobian], speye(n)}; + else + options.Jacobian = {[- options.Jacobian], eye(n)}; + endif + endif endif endif @@ -589,6 +597,16 @@ %! [t, y] = ode15s (@rob, [0, 100], [1; 0; 0], opt); %! assert ([t(end), y(end,:)], frefrob, 1e-3); +## Jacobian as const matrix +%!testif HAVE_SUNDIALS +%! opt = odeset ("RelTol", 1e-4, "AbsTol", 1e-5, +%! "Jacobian", [98, 198; -99, -199]); +%! [t, y] = ode15s (@(t, y)[98, 198; -99, -199] * (y - [1; 0]), +%! [0, 5], [2; 0], opt); +%! y1xct = @(t) 2 * exp (-t) - exp (-100 * t) + 1; +%! y2xct = @(t) - exp (-t) + exp (-100 * t); +%! assert ([y1xct(t), y2xct(t)], y, 1e-3); + ## two output arguments %!testif HAVE_SUNDIALS %! [t, y] = ode15s (@fpol, [0, 2], [2, 0]);