find numerical solution of stiff ordinary differential equations
dsolve(odesys, numeric, method=rosenbrock, vars, options)
dsolve(numeric, method=rosenbrock, procopts, options)
set or list; ordinary differential equation(s) and initial conditions
literal name; instruct dsolve to find a numerical solution
literal equation; numerical method to use
(optional) any indeterminate function of one variable, or a set or list of them, representing the unknowns of the ODE problem
(optional) equations of the form keyword = value
options used to specify the ODE system using a procedure (procedure, initial, start, number, procvars, and jacobian). For more information on all but the last of these, see dsolve/numeric/IVP.
The dsolve command with options numeric and method=rosenbrock finds a numerical solution using an Implicit Rosenbrock third-fourth order Runge-Kutta method with degree three interpolant. This is the default method of the type=numeric solution for initial value problems when the stiff argument is used.
Modes of Operation
The rosenbrock method has two distinct modes of operation (for procedure-type outputs).
With the range option specified
When used with the range option, the method computes the solution for the IVP over the specified range, storing the solution information internally, and uses that information to rapidly interpolate the desired solution value for any call to the returned procedure(s).
Though possible, it is not recommended that the returned procedure be called for points outside the specified range.
This method can be used in combination with the refine option of odeplot to produce an adaptive plot (that is, odeplot uses the precomputed points to produce the plot when refine is specified.)
It is not recommended that this method be used for problems in which the solution can become singular, as each step is stored, and many steps may be taken when near a singularity, so memory usage can become a significant issue.
The storage of the interpolant in use by this method can be disabled by using the interpolation=false option described below. This may be desirable for high accuracy solutions where storage of the interpolant (in addition to the discrete solution) requires too much memory. Disabling the interpolant is not generally recommended as the solution values will be obtained from an interpolation of the 4 closest points, and does not necessarily provide an interpolant with order 3 error.
Without the range option specified
When used without the range option, the IVP solution values are not stored, but rather computed when requested.
Because not all solution values are stored, computation must restart at the initial values whenever a point is requested between the initial point and the most recently computed point (to avoid reversal of the integration direction), so it is advisable to collect solution values moving away from the initial value.
The following options are available for the rosenbrock method.
keyword or array
name or list of names
boolean or auto
Specifies the desired output from dsolve. The keywords procedurelist, listprocedure, or operator provide procedure-type output, the keyword piecewise provides output in the form of piecewise functions over a specified range of independent variable values, and a 1-D array or Array provide output at fixed values of the independent variable. For details, see dsolve/numeric.
Specifies user-defined known functions, and basic usage is discussed in dsolve/numeric. As mentioned there, in addition to the requirement of a procedural definition for the function that evaluates with numeric arguments, otherwise not, this method also requires a `diff/` rule.
This `diff/` rule must provide the partial derivatives of the function with respect to the independent variables, dependent variables, and derivatives that are used to compute the known function. Note: The key word here is partial, and an example can be found below. This is in contrast to the taylorseries method, for which the diff rule with respect to the independent variable must provide a total derivative.
Only a single partial derivative with respect to the independent variables, dependent variables, and derivatives is required for each known function, so these can be coded as procedures themselves, not requiring `diff/` rules.
abserr, relerr, and initstep
Specify the desired accuracy of the solution, and the starting step size for the method, and are discussed in dsolve/Error_Control. The default value for rosenbrock for both the abserr and relerr options is 1e-6. The value for initstep, if not specified, is determined by the method, taking into account the local behavior of the ODE system.
By default this is set to true, and controls whether the solution interpolant error (including the interpolant on index-1 variables for DAE problems) is integrated into the error control. When set to false, areas where the solutions is varying rapidly (e.g. a discontinuity in a derivative due to a piecewise) may have a much larger solution error than dictated by the specified error tolerances. When set to true, the step size is reduced to minimize error in these regions, but for problems where there is a jump discontinuity in the variables, the integration may fail with an error indicating that a singularity is present. In the latter case where an error is thrown, it may be advantageous to model the discontinuities using events (see dsolve/Events).
Specifies the maximum number of evaluations of the right-hand side of the first order ODE system. This option can be disabled by specifying maxfun=0. The default value for rosenbrock is 30000.
number, procedure, start, initial, and procvars
These options are used to specify the IVP using procedures. For more information, see dsolve/numeric/IVP.
This option can be used in combination with the procedure and related options to specify a procedure-defined jacobian. If not provided, a jacobian is computed from the input system, or using the input procedure by evaluating it with indeterminate values of the independent and dependent variables, and utilizing the result. Note that this approach will not work if the input procedure-form system contains conditional statements and/or loops, so in these cases use of a specified procedure-form jacobian will be necessary.
As with the procedure option, the IVP system must first be written as a first order system (see dsolve/numeric/IVP), and from this the jacobian can be obtained. The input procedure must be of the form jac⁡x,y,fx,fy, where x is the value of the independent variable (input), y is the vector of values of the dependent variables (input), fx is a vector of values of the derivatives of the right-hand-sides of the differential system with respect to x (output), and fy is a matrix of values of the derivatives of the right-hand-sides of the differential system with respect to each dependent variable y (output).
As an example, consider the differential system ⅆⅆtx⁡t=f⁡t,x⁡t,y⁡t,ⅆⅆty⁡t=g⁡t,x⁡t,y⁡t for some specific form of f⁡t,x⁡t,y⁡t,g⁡t,x⁡t,y⁡t. The procedure and jacobian procedures could be defined as follows:
prc := proc(n,t,v,vp)
vp := f(t,v,v);
vp := g(t,v,v);
jac := proc(t,v,ft,fv)
ft := D(f)(t,v,v); ft := D(g)(t,v,v);
fv[1,1] := D(f)(t,v,v); fv[1,2] := D(f)(t,v,v);
fv[2,1] := D(g)(t,v,v); fv[2,2] := D(g)(t,v,v);
startinit,implicit, and optimize
These options control the method and behavior of the computation. For more information on the first two, see dsolve/numeric/IVP, for the last, see dsolve/numeric.
This option specifies that the internally generated procedures that are used to compute the numeric solution be compiled for efficiency. Note that this option will only work if Digits is set within the hardware precision range, and the input function contains only evalhf capable functions (e.g. only elementary mathematical functions like exp, sin, and ln). See dsolve/Efficiency. By default this value is set to false. If set to true and a compile is not possible, an error will be thrown. If set to auto and a compile is not possible, the uncompiled procedures will be used directly.
Determines the range of values of the independent variable for which solution values are required. Use of this option significantly changes the behavior of the method for the procedure-style output types discussed in dsolve/numeric (see description of with/without range above).
events, event_pre, event_maxiter, event_iterate, event_initial
These options are used to specify and control the behavior of events for the numerical solution. These options are discussed in detail in dsolve/Events.
Accepts a boolean value indicating if the problem is (or will become) complex valued. By default this is detected based on the input system and initial data, but in cases where the input system is procedure defined, or the system is initially real, it may be necessary to specify complex=true to obtain the solution. It is assumed that for an initially real system that becomes complex, the point at which this transition occurs is considered to be a singularity, so if complex=true is not specified, the integration will halt at that point.
By setting infolevel[dsolve] to 2, information on the last mesh point evaluation is provided in the event of an error.
The rosenbrock method is capable of computing high-accuracy solutions for IVPs because the precision of the computation can be increased by changing the Digits environment variable. Care must be taken because computation of high accuracy solutions with a method that has not been designed for the task can be quite inefficient.
Results can be plotted by using the function odeplot in the plots package.
A stiff linear problem with eigenvalues of −1, −100, storing solution values from 0..10:
deq1 ≔ ⅆ2ⅆt2⁢y⁡t+101⁢ⅆⅆt⁢y⁡t+100⁢y⁡t=0:
ic1 ≔ y⁡0=1,D⁡y⁡0=−1:
dsol1 ≔ dsolve⁡deq1,ic1,numeric,stiff=true,range=0..10
dsol1 ≔ procx_rosenbrock...end proc
A stiff nonlinear system from reactor kinetics, computing solution values as requested:
deq2 ≔ ⅆⅆt⁢u⁡t=0.01−1+u⁡t+1000⁢u⁡t+1⁢0.01+u⁡t+v⁡t,ⅆⅆt⁢v⁡t=0.01−1+v⁡t2⁢0.01+u⁡t+v⁡t
ic2 ≔ u⁡0=0,v⁡0=0:
dsol2 ≔ dsolve⁡deq2,ic2,numeric,method=rosenbrock
dsol2 ≔ procx_rosenbrock...end proc
Use of the known option for multiple argument functions (requires diff rule). Using the following function
define 'f' as
f := proc(x,y) local t;
if not type(evalf(x),numeric) or
not type(evalf(y),numeric) then
f ≔ procx,ylocalt;ifnottype⁡evalf⁡x,numericandtype⁡evalf⁡y,numericthen'procname'⁡x,yelseevalf⁡Int⁡1/10*exp⁡−1/50*t,t=0..x − yend ifend proc
and use the `diff/` rule for the partial derivatives with respect to both x and y:
`diff/f` := proc(x,y);
if args=x then
elif args=y then
error "unable to differentiate";
diff/f ≔ procx,yifargs=xthen1/10*diff⁡x,args*exp⁡−1/50*xelifargs=ythen−1elseerrorunable to differentiateend ifend proc
Using this information you can now call dsolve to obtain a solution procedure, then obtain solution values.
dsol ≔ dsolve⁡ⅆⅆx⁢y⁡x=f⁡x,y⁡x,y⁡0=0,numeric,known=f,method=rosenbrock
dsol ≔ procx_rosenbrock...end proc
Hairer, E., and Wanner, G. Solving Ordinary Differential Equations II. 2nd ed. New York: Springer, 1996.
Shampine, L.F., and Corless, R.M. "Initial Value Problems for ODEs in Problem Solving Environments." J. Comp. Appl. Math, Vol. 125(1-2), (2000): 31-40.
Download Help Document