find numerical solution of ODE initial value problems - Maple Help

Online Help

All Products    Maple    MapleSim


Home : Support : Online Help : Mathematics : Differential Equations : dsolve : numeric : dsolve/numeric/IVP

dsolve/numeric/IVP - find numerical solution of ODE initial value problems

Calling Sequence

dsolve(odesys, numeric, vars, options)

dsolve(numeric, procopts, options)

Parameters

odesys

-

set or list; ordinary differential equation(s) and initial conditions

numeric

-

name; instruct dsolve to find a numerical solution

vars

-

(optional) any indeterminate function of one variable, or a set or list of them, representing the unknowns of the ODE problem

procopts

-

(required if odesys is not present) options that specify a procedure-defined system (procedure, initial, start, number, and procvars). For more information, see below.

options

-

(optional) equations of the form keyword = value

Description

• 

The dsolve command with the numeric or type=numeric option and an initial value problem (IVP) finds a numerical solution for the ODE or ODE system IVP. If the optional equation method=numericmethod is provided (where numericmethod is one of rkf45, ck45, rosenbrock, dverk78, lsode, gear, taylorseries, or classical), dsolve uses that method to obtain the numerical solution.

• 

The return value of dsolve and the following high level or common options are discussed in dsolve[numeric] and dsolve[Error_Control].

'output'     =

keyword or array

'range'      =

numeric..numeric

'abserr'     =

numeric or list

'relerr'     =

numeric

'stiff'      =

boolean

'known'      =

name or list of names

• 

The exception is the list input form for abserr which requires the discussion below on conversion to first order to understand fully.

  

In most cases it is sufficient to specify a single value for abserr for all solution components. The exception to this is when some dependent variables for derivatives of the system have values on a different scale than other variables of the system. When this occurs, the list input form for abserr can be used to specify a different absolute error tolerance for each dependent variable of the system, or each solution component of the system (the difference between these is that specification of a tolerance for a dependent variable of the system also does so for derivatives of that variable, while for solution components each value is considered independently.

  

In order to specify dependent or component abserr, the order of the dependent variables in the output must be known. This can be controlled by specifying the dependent variables of the problem (the vars parameter to dsolve) as a list.

  

From this, dependent variable abserr is easy to understand as there is a 1-1 correspondence between the specified variables and the specified abserr.

  

The per-component abserr is a little trickier as the basic idea is the same, but depending on the order of each variable more than one abserr value may be required.

  

For example, for a problem containing x'',y''',z', if the variables were specified as xt,yt,zt, then for per-component abserr six values would be required, and would correspond to [xt,x't,yt,y't,y''t,zt].

• 

Only the rkf45, ck45, and rosenbrock IVP methods support the list form for abserr.

• 

The default IVP method is a Runge-Kutta Fehlberg method, which produces a solution accurate to fifth order. Another non-stiff IVP method is a Runge-Kutta method with the Cash-Karp coefficients. The default stiff method is a Rosenbrock method, which produces a solution accurate to fourth order. See dsolve[rkf45], dsolve[ck45], and dsolve[rosenbrock] for more detail on these methods. The other available methods for IVP are dsolve[dverk78], dsolve[lsode], dsolve[taylorseries], dsolve[gear], and dsolve[classical].

• 

All IVP methods, with the exception of taylorseries, can work with a complex-valued IVP with a real-valued independent variable.

• 

The following options are specific to IVP, and are concerned with specification of the input system in procedure form (procopts):

'number'    =

integer

'procedure' =

procedure

'start'     =

numeric

'initial'   =

array

'procvars'  =

list

  

and the following options, also specific to IVP (and in some cases DAE), allow control of the solution method and step size, and the use of parameters in the system:

'startinit'     =

boolean

'minstep'       =

numeric

'maxstep'       =

numeric

'initstep'      =

numeric

'events'        =

list

'event_pre'     =

keyword

'event_maxiter' =

integer

'implicit'      =

boolean

'parameters'    =

list of names

'complex'       =

boolean

• 

Before discussion of the procedure-based options, the Maple treatment of the numerical handling of ODE system IVPs is described. To begin with, all systems are converted to first order systems by the introduction of new dependent variables. For example, if the single ODE

3y'''=2y'xy=0

  

were to be numerically solved, the new dependent variables u=y', v=y'' would be introduced, and the single third order ODE would become a system of three first order ODEs in three dependent variables. In addition, the derivatives of the first order system are isolated. The result would be

y'=u

u'=v

v'=13xy+23u

  

and the initial conditions would transform in the same way.

  

Now that the ODE is converted to a first order system, it is possible to describe the procedure-based options more precisely. In all that follows, any reference to the system or dependent variables refers only to the first order system.

Options

  

'number'= integer

  

The number of dependent variables in the first order system.

  

 

  

'procedure'= procedure

  

Procedure for the evaluation of the derivatives of the dependent variables. It must be of the form procN,x,Y,YP, where N is the number of dependent variables, x is the value of the independent variable at which the derivatives are to be evaluated, Y is an array of the dependent variable values at x, and YP is an output array into which the procedure places the computed values of the derivatives.

  

For the above example, given that the correspondence of the variables is Y1=y,Y2=u, and Y3=v, the following procedure correctly evaluates the derivatives of the ODE system.

solproc := proc(N, x, Y, YP)

   YP[1] := Y[2];

   YP[2] := Y[3];

   YP[3] := (2*Y[2] - x*Y[1])/3;

end proc;

  

 

  

'start'= numeric

  

Numeric value that specifies the initial value of the independent variable.

  

 

  

'initial'= array

  

Array that specifies the initial values of the dependent variables of the first order system. These values must be present in this array in the same order as that used by the procedure given in the procedure argument.

  

 

  

'procvars'= list

  

List that specifies the names of the dependent variables corresponding to the indexed values in the solution procedure and the initial value array. For the above example, procvars=yx,ux,vx can be used. To instead display solutions in the variables of the original ODE, procvars=yx,ⅆⅆxyx,ⅆ2ⅆx2yx can be used.

  

When an ODE or ODE system is input directly (a system defined problem - does not use the procedure options), Maple converts the system to first order internally by renaming derivatives as new dependent variables. It obtains the number, procedure, start, initial, and procvars values from the input. In cases where both the system form and the procedure form of the input is present, the system form is ignored and a warning is issued.

  

Note: If the ODE system is provided in solved form with respect to the leading derivatives (the highest derivative of each dependent variable of the system), no algebraic manipulation of the system is done in obtaining the procedure form (other than derivative renaming to new dependent variables). This allows better control of the actual computation of the derivatives in the system. If any of the equations in the ODE system are not in solved form, solve is called to obtain the solved form, which means that certain undesirable side effects may occur (for example, consider solve(y-(x+1)^20-1,y) ).

  

 

  

'implicit'= boolean

  

Argument that can only be true for ODE systems that are linear in the leading derivatives, and can only be used with the rkf45, ck45, rosenbrock, dverk78, gear, lsode, and classical methods with system-defined problems. By default, this option is false, but setting this option to true tells dsolve[numeric] to leave the system in unsolved form, and compute the derivative values through a linear solve each time they are needed. This means that a method requiring n function evaluations per step must perform n linear solves per step.

  

Though use of this option slows the computation in most cases, it is necessary for ODE systems that have complicated coefficients, for which the isolation of the leading derivatives would result in excessive expression swell, or an ill-conditioned function to compute the derivative values. This behavior occurs frequently in mapping solutions of polynomial systems to solutions of differential systems using homotopy methods.

  

 

  

'events'= list

  

'event_pre'= keyword

  

'event_maxiter'= integer

  

Arguments that are only available for the rkf45, ck45, and rosenbrock methods, and are used to specify and control the behavior of events. These options are discussed in detail in dsolve[Events].

  

 

  

'startinit'= boolean

  

Boolean that controls the behavior of the numerical integration with respect to consecutive calculations. If startinit=true is specified, each calculation begins at the initial point of the problem. If startinit=false (the default setting), then each calculation begins where the prior calculation ended, provided it does not reverse the direction of the integration. If doing so would reverse the integration direction, the computation is started from the initial point.

  

For the methods, rkf45, ck45, and rosenbrock, the solution is only recomputed when the new value of the independent variable is not between the initial point and any other previously requested solution point. This has the effect of never reversing the direction of integration, and making evaluation of the solution for an already computed interval quite inexpensive (the solution values are obtained by interpolation). The storage of the solution can be disabled by using the storage argument. The startinit parameter also tells these methods to recompute the solution whenever a solution value is desired.

  

 

  

'complex'=boolean

  

For all IVP methods, you may specify that the problem is complex with the complex option.  Note that this option is not available for the other numeric methods. For the default stiff and nonstiff IVP methods, it may be necessary to specify when the problem is complex.  The default value is false.

  

 

• 

The minstep, maxstep, and initstep options are discussed in dsolve[Error_Control].

• 

In addition to computation of solution values for the given input problem, the procedure returns (that is, output=procedurelist, listprocedure or operator) provide additional interactive features, including the ability to specify parameters. Information on these features is provided on the dsolve[numeric,interactive] page.

Examples

Initial value problem with solution yx=sinx, integrated in both directions:

dsys1:=ⅆ2ⅆx2yx+yx=0,y1=sin1,Dy1=cos1

dsys1:=ⅆ2ⅆx2yx+yx=0,y1=sin1,Dy1=cos1

(1)

dsol1:=dsolvedsys1,numeric,method=gear

dsol1:=procx_gear...end proc

(2)

dsol10

x=0.,yx=8.0533900612613710-12,ⅆⅆxyx=0.999999999997424

(3)

dsol1π

x=3.14159265358979,yx=6.9771000387735610-12,ⅆⅆxyx=0.999999999999881

(4)

IVP system:

deq2:=ⅆⅆtxt=yt,ⅆⅆtyt=xt+yt

deq2:=ⅆⅆtxt=yt,ⅆⅆtyt=xt+yt

(5)

ic2:=x0=2,y0=1:

dsol2:=dsolvedeq2∪ic2,numeric,output=listprocedure,range=0..1

dsol2:=t=proct...end proc,xt=proct...end proc,yt=proct...end proc

(6)

dsol2x:=subsdsol2,xt:dsol2y:=subsdsol2,yt:

dsol2x0,dsol2y0

2.,1.

(7)

dsol2x0.4,dsol2y0.4

2.69118458493332,2.60811748317261

(8)

dsol2x1.0,dsol2y1.0

5.58216753140384,7.82688926499912

(9)

Lorenz system using procedure specification:

dproc3 := proc(N,t,Y,YP)
   YP[1] := sigma * (Y[2] - Y[1]);
   YP[2] := r*Y[1] - Y[1]*Y[3] - Y[2];
   YP[3] := Y[1]*Y[2] - b*Y[3];
end proc:

σ:=10:r:=28:b:=83:

ic3:=Array5,5,5:

dvars3:=xt,yt,zt:

dsol3:=dsolvenumeric,number=3,procedure=dproc3,start=0,initial=ic3,procvars=dvars3

dsol3:=procx_rkf45...end proc

(10)

dsol31

t=1.,xt=7.09064281000129,yt=4.13869080943593,zt=29.0616065518357

(11)

dsol32

t=2.,xt=9.13349830723850,yt=12.6955376022149,zt=22.4637998042670

(12)

dsol33

t=3.,xt=5.00807269946467,yt=2.53517706788747,zt=26.6155378375388

(13)

The following is a problem (of the form ⅇxy''+sinxyx=0 ) where the natural naming of the variables is not the desired one. Choose the first dependent variable as yx, and the second as vx=ⅇxy'x, and use the procedure specification of the problem.

dproc4 := proc(N,x,Y,YP)
   YP[1] := exp(-x)*Y[2];
   YP[2] := -sin(x)*Y[1];
end proc;

dproc4:=procN,x,Y,YPYP[1]:=exp−x*Y[2];YP[2]:=−sinx*Y[1]end proc

(14)

ic4:=Array1,0:

dvars4:=yx,ⅇxⅆⅆxyx:

dsol4:=dsolvenumeric,number=2,procedure=dproc4,start=0,initial=ic4,procvars=dvars4

dsol4:=procx_rkf45...end proc

(15)

dsol40

x=0.,yx=1.,ⅇxⅆⅆxyx=0.

(16)

dsol41

x=1.,yx=0.924554126889498,ⅇxⅆⅆxyx=0.444194433465032

(17)

dsol42

x=2.,yx=0.741972301895190,ⅇxⅆⅆxyx=1.24013689804566

(18)

A complex-valued initial value problem:

dsys1:=ⅆⅆxyx+Iyx=0,y0=1+2I

dsys1:=ⅆⅆxyx+Iyx=0,y0=1+2I

(19)

dsol5:=dsolvedsys1,numeric,method=gear

dsol5:=procx_gear...end proc

(20)

dsol5π2

x=1.57079632679490,yx=2.000000000003670.999999999983486I

(21)

dsol5π

x=3.14159265358979,yx=1.000000000014092.00000000003328I

(22)

dsol52π

x=6.28318530717958,yx=1.00000000000419+2.00000000002333I

(23)

Interactive examples:

dsys:=ⅆⅆxyx=yx2,y0=1

dsys:=ⅆⅆxyx=yx2,y0=1

(24)

dsoli:=dsolvedsys,numeric

dsoli:=procx_rkf45...end proc

(25)

Perform a computation.

dsoli0.5

x=0.5,yx=2.00000045906252

(26)

Query the initial value, last computed value, and method.

dsoli'initial'

x=0.,yx=1.

(27)

dsoli'last'

x=0.500000000000000,yx=2.00000045906252

(28)

dsoli'method'

rkf45

(29)

Compute up to a singularity.

dsoli1

Error, (in dsoli) cannot evaluate the solution further right of .99999999, probably a singularity

Query closest point to singularity that the method obtained.

dsoli'last'

x=0.999999994312848,yx=1.006071915623691012

(30)

Change initial condition, and recompute.

dsoli'initial'=1,12

x=1.,yx=0.500000000000000

(31)

dsoli3

Error, (in dsoli) cannot evaluate the solution further right of 2.9999999, probably a singularity

dsoli'last'

x=2.99999998871727,yx=3.257285711635441011

(32)

Now listprocedure output:

dsoli2:=dsolvedsys,numeric,output=listprocedure

dsoli2:=x=procx...end proc,yx=procx...end proc

(33)

dsx:=rhsdsoli21

dsx:=procx...end proc

(34)

dsy:=rhsdsoli22

dsy:=procx...end proc

(35)

Compute and query.

dsy0.5

2.00000045906252

(36)

dsx'last'

0.500000000000000

(37)

Change the initial x point only, and compute y1.5.

dsx'initial'=1

1.

(38)

dsy1.5

2.00000045906252

(39)

Change both.

dsy'initial'=1,12

0.500000000000000

(40)

dsy3.0

Error, (in dsy) cannot evaluate the solution further right of 2.9999999, probably a singularity

dsx'last',dsy'last'

2.99999998871727,3.257285711635441011

(41)

See Also

dsolve/Error_Control, dsolve[ck45], dsolve[classical], dsolve[dverk78], dsolve[Events], dsolve[gear], dsolve[lsode], dsolve[maxfun], dsolve[numeric,DAE], dsolve[numeric,interactive], dsolve[numeric], dsolve[rkf45], dsolve[rosenbrock], dsolve[Stiffness], dsolve[taylorseries], plots[odeplot]


Download Help Document

Was this information helpful?



Please add your Comment (Optional)
E-mail Address (Optional)
What is ? This question helps us to combat spam