Using MATLAB (779505), страница 55
Текст из файла (страница 55)
It is a one-step solver – in computing y(t n) , itneeds only the solution at the immediately preceding time point,y(t n – 1) . In general, ode45 is the best function to apply as a “firsttry” for most problems.ode23Based on an explicit Runge-Kutta (2,3) pair of Bogacki andShampine.
It may be more efficient than ode45 at crude tolerancesand in the presence of mild stiffness. Like ode45, ode23 is aone-step solver.ode113Variable order Adams-Bashforth-Moulton PECE solver. It may bemore efficient than ode45 at stringent tolerances and when theODE function is particularly expensive to evaluate. ode113 is amultistep solver – it normally needs the solutions at severalpreceding time points to compute the current solution.Solvers for Stiff ProblemsNot all difficult problems are stiff, but all stiff problems are difficult for solversnot specifically designed for them.
Solvers for stiff problems can be used exactly15-715Differential Equationslike the other solvers. However, you can often significantly improve theefficiency of these solvers by providing them with additional information aboutthe problem. (See “Improving ODE Solver Performance” on page 15-16.)There are four solvers designed for stiff problems:ode15sVariable-order solver based on the numerical differentiationformulas (NDFs).
Optionally it uses the backward differentiationformulas, BDFs, (also known as Gear’s method). Like ode113,ode15s is a multistep solver. If you suspect that a problem is stiff orif ode45 failed or was very inefficient, try ode15s.ode23sBased on a modified Rosenbrock formula of order 2. Because it is aone-step solver, it may be more efficient than ode15s at crudetolerances. It can solve some kinds of stiff problems for whichode15s is not effective.ode23tAn implementation of the trapezoidal rule using a “free”interpolant. Use this solver if the problem is only moderately stiffand you need a solution without numerical damping.ode23tb An implementation of TR-BDF2, an implicit Runge-Kutta formulawith a first stage that is a trapezoidal rule step and a second stagethat is a backward differentiation formula of order 2.
Like ode23s,this solver may be more efficient than ode15s at crude tolerances.ODE Solver Basic SyntaxAll of the ODE solver functions share a syntax that makes it easy to try any ofthe different numerical methods if it is not apparent which is the mostappropriate. To apply a different method to the same problem, simply changethe ODE solver function name. The simplest syntax, common to all the solverfunctions, is[t,y] = solver(odefun,tspan,y0)where solver is one of the ODE solver functions listed previously.15-8Initial Value Problems for ODEs and DAEsThe basic input arguments are:odefunFunction that evaluates the system of ODEs.
It has the formdydt = odefun(t,y)where t is a scalar, and dydt and y are column vectors.tspanVector specifying the interval of integration. The solver imposesthe initial conditions at tspan(1), and integrates from tspan(1) totspan(end).For tspan vectors with two elements [t0 tf], the solver returnsthe solution evaluated at every integration step.
For tspan vectorswith more than two elements, the solver returns solutionsevaluated at the given time points. The time values must be inorder, either increasing or decreasing.Specifying tspan with more than two elements does not affect theinternal time steps that the solver uses to traverse the intervalfrom tspan(1) to tspan(end). All solvers in the MATLAB ODEsuite obtain output values by means of continuous extensions ofthe basic formulas. Although a solver does not necessarily stepprecisely to a time point specified in tspan, the solutions producedat the specified time points are of the same order of accuracy as thesolutions computed at the internal time points.Specifying tspan with more than two elements has little effect onthe efficiency of computation, but for large systems, affectsmemory management.y0Vector of initial conditions for the problemSee also “Introduction to Initial Value ODE Problems” onpage 15-6.The output arguments are:tColumn vector of time pointsySolution array.
Each row in y corresponds to the solution at a timereturned in the corresponding row of t.15-915Differential EquationsAdditional ODE Solver ArgumentsFor more advanced applications, you can also specify as input arguments solveroptions and additional problem parameters.optionsStructure of optional parameters that change the defaultintegration properties. This is the fourth input argument.[t,y] = solver(odefun,tspan,y0,options)“Improving ODE Solver Performance” on page 15-16 tells youhow to create the structure and describes the properties youcan specify.p1,p2...Parameters that the solver passes to odefun.[t,y] = solver(odefun,tspan,y0,options,p1,p2...)The solver passes any input parameters that follow the optionsargument to odefun and any functions you specify in options.Use options = [] as a placeholder if you set no options.
Thefunction odefun must have the formdydt = odefun(t,y,p1,p2,...)See “Passing Additional Parameters to an ODE Function” onpage 15-13 for an example.Representing ODE ProblemsThis section describes the process for solving initial value ODE problems usingone of the MATLAB ODE solvers. It uses the van der Pol equation:• To describe the steps needed to solve an ODE problem using MATLAB• As an example of a stiff ODE problemIt also tells you how to evaluate the solution at specific points.Note See “ODE Solver Basic Syntax” on page 15-8 for more information.15-10Initial Value Problems for ODEs and DAEsExample: Solving an IVP ODE in MATLAB (van der Pol Equation, Nonstiff)This example explains and illustrates the steps you need to solve an initialvalue ODE problem using MATLAB.1 Rewrite the Problem as a System of First-Order ODEs. ODEs often involve a numberof dependent variables, as well as derivatives of order higher than one.
To usethe MATLAB ODE solvers, you must rewrite such equations as an equivalentsystem of first-order differential equations of the formy′ = f ( t, y )You can write any ordinary differential equationy(n)= f ( t, y, y′, …, y(n – 1))as a system of first-order equations by making the substitutionsy 1 = y, y 2 = y′, …, y n = y(n – 1)The result is an equivalent system of n first-order ODEs.y′1 = y 2y 2′ = y 3…y n′ = f ( t, y 1, y 2, ..., y n )For example, you can rewrite the van der Pol equation (second-order)2y′′1 – µ ( 1 – y 1 ) y′1 + y 1 = 0where µ > 0 is a scalar parameter, by making the substitution y′ 1 = y 2 .
Theresulting system of first-order ODEs isy′ 1 = y 22y′ 2 = µ ( 1 – y 1 )y 2 – y 12 Code the System of First-Order ODEs in MATLAB. Once you represent the equationas a system of first-order ODEs, you can code it as a function that a MATLABODE solver can use. The function must be of the form15-1115Differential Equationsdydt = odefun(t,y)Although t and y must be the function’s first two arguments, the function doesnot need to use them. The output dydt, a column vector, is the derivative of y.The code below represents the van der Pol system in a MATLAB function, vdp1.The vdp1 function assumes that µ = 1 .
y 1 and y 2 become elements y(1) andy(2) of a two-element vector.function dydt = vdp1(t,y)dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];Note that, although vdp1 must accept the arguments t and y, it does not use tin its computations.3 Apply a Solver to the Problem. Decide which solver you want to use to solve theproblem. Then call the solver and pass it the function you created to describethe first-order system of ODEs, the time interval on which you want to solvethe problem, and an initial condition vector. See “Initial Value ProblemSolvers” on page 15-6 and the ODE solver reference page for descriptions of theODE solvers.For the van der Pol system, you can use ode45 on time interval [0 20] withinitial values y(1) = 2 and y(2) = 0.[t,y] = ode45(@vdp1,[0 20],[2; 0]);This example uses @ to pass vdp1 as a function handle to ode45. The resultingoutput is a column vector of time points t and a solution array y.
Each row iny corresponds to a time returned in the corresponding row of t. The first columnof y corresponds to y 1 , and the second column to y 2 .Note See the function_handle (@), func2str, and str2func reference pages,and the Function Handles chapter of “Programming and Data Types” in theMATLAB documentation for information about function handles.4 View the Solver Output. You can simply use the plot command to view thesolver output.plot(t,y(:,1),'-',t,y(:,2),'--')title('Solution of van der Pol Equation, \mu = 1');15-12Initial Value Problems for ODEs and DAEsxlabel('time t');ylabel('solution y');legend('y_1','y_2')Solution of van der Pol Equation, µ = 13y1y22solution y10−1−2−30246810time t1214161820As an alternative, you can use a solver output function to process the output.The solver calls the function specified in the integration property OutputFcnafter each successful time step. Use odeset to set OutputFcn to the desiredfunction.
See “OutputFcn” on page 15-20 for more information.Passing Additional Parameters to an ODE Function. The solver passes any inputparameters that follow the options argument to the ODE function and anyfunction you specify in options. For example:1 Generalize the van der Pol function by passing it a mu parameter, instead ofspecifying a value for mu explicitly in the code.function dydt = vdp1(t,y,mu)dydt = [y(2); mu*(1-y(1)^2)*y(2)-y(1)];15-1315Differential Equations2 Pass the parameter mu to the function vdp1 by specifying it after the optionsargument in the call to the solver.
This example uses options = [] as aplaceholder.[t,y] = ode45(@vdp1,tspan,y0,[],mu)callsvdp1(t,y,mu)See vdpode for a complete example based on these functions.Example: The van der Pol Equation, µ = 1000 (Stiff)This example presents a stiff problem. For a stiff problem, solutions can changeon a time scale that is very short compared to the interval of integration, butthe solution of interest changes on a much longer time scale. Methods notdesigned for stiff problems are ineffective on intervals where the solutionchanges slowly because they use time steps small enough to resolve the fastestpossible change.When µ is increased to 1000, the solution to the van der Pol equation changesdramatically and exhibits oscillation on a much longer time scale.Approximating the solution of the initial value problem becomes a moredifficult task.















