Using MATLAB (779505), страница 56
Текст из файла (страница 56)
Because this particular problem is stiff, a solver intended fornonstiff problems, such as ode45, is too inefficient to be practical. A solver suchas ode15s is intended for such stiff problems.The vdp1000 function evaluates the van der Pol system with µ = 1000.function dydt = vdp1000(t,y)dydt = [y(2); 1000*(1-y(1)^2)*y(2)-y(1)];Note This example hardcodes µ in the ODE function. The vdpode examplesolves the same problem, but passes a user-specified µ as an additionalargument to the ODE function.
See “Additional ODE Solver Arguments” onpage 15-10.15-14Initial Value Problems for ODEs and DAEsNow use the ode15s function to solve the problem with the initial conditionvector of [2; 0], but a time interval of [0 3000]. For scaling purposes, plot justthe first component of y(t).[t,y] = ode15s(@vdp1000,[0 3000],[2; 0]);plot(t,y(:,1),'-');title('Solution of van der Pol Equation, \mu = 1000');xlabel('time t');ylabel('solution y_1');Solution of van der Pol Equation, µ = 10002.521.51solution y10.50−0.5−1−1.5−2−2.5050010001500time t200025003000Note For detailed instructions for solving an initial value ODE problem inMATLAB, see “Example: Solving an IVP ODE in MATLAB (van der PolEquation, Nonstiff)” on page 15-11.15-1515Differential EquationsEvaluating the Solution at Specific PointsThe numerical methods implemented in the ODE solvers produce a continuoussolution over the interval of integration [ a, b ] .
You can evaluate theapproximate solution, S ( x ) , at any point in [ a, b ] using the function deval andthe structure sol returned by the solver.Sxint = deval(sol,xint)The ODE solvers return the structure sol when called with a single outputargument.The deval function is vectorized. For a vector xint, the ith column of Sxintapproximates the solution y ( xint(i) ) .Improving ODE Solver PerformanceThe default integration properties in the ODE solvers are selected to handlecommon problems.
In some cases, you can improve ODE solver performance byoverriding these defaults. You do this by supplying the solvers with one or moreproperty values in an options structure.[t,y] = solver(odefun,tspan,y0,options)This section:• Explains how to create, modify, and query an options structure.• Describes the properties that you can use in an options structure.In this and subsequent property tables, the most commonly used properties arelisted first, followed by more advanced properties.ODE Property Categories15-16Properties CategoryProperty NameError controlRelTol, AbsTol, NormControlSolver outputOutputFcn, OutputSel, Refine, StatsJacobian matrixJacobian, JPattern, VectorizedStep-sizeInitialStep, MaxStepInitial Value Problems for ODEs and DAEsODE Property CategoriesProperties CategoryProperty NameMass matrix andDAEsMass, MStateDependence, MvPattern,MassSingular, InitialSlopeEvent locationEventsode15s-specificMaxOrder, BDFCreating and Maintaining an ODE Options StructureCreating an Options Structure.
The odeset function creates an options structurethat you can pass as an argument to any of the ODE solvers. To create anoptions structure, odeset accepts property name/property value pairs usingthe syntaxoptions = odeset('name1',value1,'name2',value2,...)In the resulting structure, options, the named properties have the specifiedvalues.
Any unspecified properties contain default values in the solvers. For allproperties, it is sufficient to type only the leading characters that uniquelyidentify the property name. odeset ignores case for property names.With no input arguments, odeset displays all property names and theirpossible values. It indicates defaults with {}.Modifying an Existing Options Structure. To modify an existing options structure,oldopts, useoptions = odeset(oldopts,'name1',value1,...)This sets options equal to the existing structure oldopts, overwrites anyvalues in oldopts that are respecified using name/value pairs, and adds anynew pairs to the structure.
The modified structure is returned as an outputargument. In the same way, the commandoptions = odeset(oldopts,newopts)combines the structures oldopts and newopts. In the output argument, anyvalues in the second argument overwrite those in the first argument.15-1715Differential EquationsQuerying Options. The odeget function extracts property values from an optionsstructure created with odeset.o = odeget(options,'name')This functions returns the value of the specified property, or an empty matrix[], if the property value is unspecified in the options structure.As with odeset, it is sufficient to type only the leading characters that uniquelyidentify the property name.
Case is ignored for property names.Error Control PropertiesAt each step, the solver estimates the local error e in the ith component of thesolution. This error must be less than or equal to the acceptable error, which isa function of the specified relative tolerance, RelTol, and the specified absolutetolerance, AbsTol.|e(i)| ≤ max(RelTol*abs(y(i)),AbsTol(i))For routine problems, the ODE solvers deliver accuracy roughly equivalent tothe accuracy you request.
They deliver less accuracy for problems integratedover “long” intervals and problems that are moderately unstable. Difficultproblems may require tighter tolerances than the default values. For relativeaccuracy, adjust RelTol. For the absolute error tolerance, the scaling of thesolution components is important: if |y| is somewhat smaller than AbsTol, thesolver is not constrained to obtain any correct digits in y.
You might have tosolve a problem more than once to discover the scale of solution components.Roughly speaking, this means that you want RelTol correct digits in allsolution components except those smaller than thresholds AbsTol(i). Even ifyou are not interested in a component y(i) when it is small, you may have tospecify AbsTol(i) small enough to get some correct digits in y(i) so that youcan accurately compute more interesting components15-18Initial Value Problems for ODEs and DAEsThe following table describes the error control properties. Use odeset to set theproperties.ODE Error Control PropertiesPropertyValueDescriptionRelTolPositivescalar {1e-3}A relative error tolerance that applies to all components of thesolution vector y.
It is a measure of the error relative to the sizeof each solution component. Roughly, it controls the number ofcorrect digits in all solution components except those smallerthan thresholds AbsTol(i).The default, 1e-3, corresponds to 0.1% accuracy.AbsTolPositivescalar orvector{1e-6}Absolute error tolerances that apply to the individualcomponents of the solution vector.
AbsTol(i) is a thresholdbelow which the value of the ith solution component isunimportant. The absolute error tolerances determine theaccuracy when the solution approaches zero. Even if you are notinterested in a component y(i) when it is small, you may haveto specify AbsTol(i) small enough to get some correct digits iny(i) so that you can accurately compute more interestingcomponents.If AbsTol is a vector, the length of AbsTol must be the same asthe length of the solution vector y.
If AbsTol is a scalar, the valueapplies to all components of y.NormControlon | {off}Control error relative to norm of solution. Set this property on torequest that the solvers control the error in each integrationstep with norm(e) <= max(RelTol*norm(y),AbsTol). Bydefault the solvers use a more stringent component-wise errorcontrol.15-1915Differential EquationsSolver Output PropertiesThe solver output properties let you control the output that the solversgenerate.
Use odeset to set these properties.ODE Solver Output PropertiesPropertyValueDescriptionOutputFcnFunction{odeplot}Installable output function. The solver calls this function afterevery successful integration step.For example,options = odeset('OutputFcn',@myfun)sets the OutputFcn property to an output function, myfun, that canbe passed to an ODE solver.The output function must be of the formstatus = myfun(t,y,flag)The solver calls the specified output function with the followingflags. Note that the syntax of the call differs with the flag. Thefunction must respond appropriately:init15-20The solver calls myfun(tspan,y0,'init') beforebeginning the integration to allow the output functionto initialize.
tspan and y0 are the input arguments tothe ODE solver.Initial Value Problems for ODEs and DAEsODE Solver Output Properties (Continued)PropertyValueDescription{none}The solver calls status = myfun(t,y) after eachintegration step on which output is requested. tcontains points where output was generated during thestep, and y is the numerical solution at the points in t.If t is a vector, the ith column of y corresponds to theith element of t.When length(tspan) > 2 the output is produced atevery point in tspan.
When length(tspan) = 2 theoutput is produced according to the Refine option.myfun must return a status output value of 0 or 1. Ifstatus = 1, the solver halts integration. You can usethis mechanism, for instance, to implement a Stopbutton.doneThe solver calls myfun([],[],'done') whenintegration is complete to allow the output function toperform any cleanup chores.You can use these general purpose output functions or you can editthem to create your own. Type help function at the command linefor more information.• odeplot – time series plotting (default when you call the solverwith no output arguments and you have not specified an outputfunction)• odephas2 – two-dimensional phase plane plotting• odephas3 – three-dimensional phase plane plotting• odeprint – print solution as it is computedNote If you call the solver with no output arguments, the solverdoes not allocate storage to hold the entire solution history.15-2115Differential EquationsODE Solver Output Properties (Continued)PropertyValueDescriptionOutputSelVector ofindicesVector of indices specifying which components of the solutionvector are to be passed to the output function.
















