Using MATLAB (779505), страница 60
Текст из файла (страница 60)
The Robertson problem coded in hb1ode is a classic testproblem for codes that solve stiff ODEs.15-4715Differential Equations4y′1 = – 0.04 y 1 + 10 y 2 y 37 24y′2 = 0.04y 1 – 10 y 2 y 3 – 3 ⋅ 10 y 27 2y′3 = 3 ⋅ 10 y 2Note The Robertson problem appears as an example in the prolog toLSODI [4].In hb1ode, the problem is solved with initial conditions y 1(0) = 1 , y 2(0) = 0 ,y 3(0) = 0 to steady state. These differential equations satisfy a linearconservation law that is used to reformulate the problem as the DAE4y′1 = – 0.04 y 1 + 10 y 2 y 347 2y′2 = 0.04y 1 – 10 y 2 y 3 – 3 ⋅ 10 y 20 = y1 + y2 + y3 – 1Obviously these equations do not have a solution for y ( 0 ) with componentsthat do not sum to 1.
The problem has the form of My′ = f ( t, y ) withM =1 0 00 1 00 0 0M is obviously singular, but hb1dae does not inform the solver of this. Thesolver must recognize that the problem is a DAE, not an ODE. Similarly,although consistent initial conditions are obvious, the example uses an–3inconsistent value y 3(0) = 10 to illustrate computation of consistent initialconditions.To run this example from the MATLAB Help browser, click on the examplename. Otherwise, type hb1dae at the command line.
Note that hb1dae:15-48Initial Value Problems for ODEs and DAEs• Imposes a much smaller absolute error tolerance on y 2 than on the othercomponents. This is because y 2 is much smaller than the other componentsand its major change takes place in a relatively short time.• Specifies additional points at which the solution is computed to more clearlyshow the behavior of y 2 .• Multiplies y 2 by 104 to make y 2 visible when plotting it with the rest of thesolution.• Uses a logarithmic scale to plot the solution on the long time interval.function hb1dae%HB1DAE Stiff differential-algebraic equation (DAE)% A constant, singular mass matrixM = [1 0 00 1 00 0 0];% Use an inconsistent initial condition to test initialization.y0 = [1; 0; 1e-3];tspan = [0 4*logspace(-6,6)];% Use the LSODI example tolerances.
The 'MassSingular' property% is left at its default 'maybe' to test the automatic detection% of a DAE.options = odeset('Mass',M,'RelTol',1e-4,...'AbsTol',[1e-6 1e-10 1e-6],’Vectorized','on');[t,y] = ode15s(@f,tspan,y0,options);y(:,2) = 1e4*y(:,2);semilogx(t,y);ylabel('1e4 * y(:,2)');title(['Robertson DAE problem with a Conservation Law, '...'solved by ODE15S']);xlabel('This is equivalent to the stiff ODEs coded in HB1ODE.');% --------------------------------------------------------------15-4915Differential Equationsfunction out = f(t,y)out = [ -0.04*y(1,:) + 1e4*y(2,:).*y(3,:)0.04*y(1,:) - 1e4*y(2,:).*y(3,:) - 3e7*y(2,:).^2y(1,:) + y(2,:) + y(3,:) - 1 ];Robertson DAE problem with a Conservation Law, solved by ODE15S10.90.80.71e4 * y(:,2)0.60.50.40.30.20.10−610−410−202410101010This is equivalent to the stiff ODEs coded in HB1ODE.610810Questions and Answers, and TroubleshootingThis section contains a number of tables that answer questions about the useand operation of the MATLAB ODE solvers:• General ODE solver questions• Problem size, memory use, and computation speed• Time steps for integration• Error tolerance and other options• Solving different kinds of problems• Troubleshooting15-50Initial Value Problems for ODEs and DAEsGeneral ODE Solver QuestionsQuestionAnswerHow do the ODE solversdiffer from quad or quadl?quad and quadl solve problems of the form y′ = f ( t ) .
The ODEsolvers handle more general problems y′ = f ( t, y ) , or problemsthat involve a mass matrix M(t, y) y′ = f(t, y) .Can I solve ODE systems inwhich there are moreequations than unknowns,or vice versa?No.Problem Size, Memory Use, and Computation SpeedQuestionAnswerHow large a problem can Isolve with the ODE suite?The primary constraints are memory and time. At each time step,the solvers for nonstiff problems allocate vectors of length n,where n is the number of equations in the system.
The solvers forstiff problems allocate vectors of length n, but also an n-by-nJacobian matrix. For these solvers it may be advantageous to usethe sparse option.If the problem is nonstiff, or if you are using the sparse option, itmay be possible to solve a problem with thousands of unknowns.In this case, however, storage of the result can be problematic.Try asking the solver to evaluate the solution at specific pointsonly, or call the solver with no output arguments and use anoutput function to monitor the solution.I’m solving a very largesystem, but only care abouta couple of the componentsof y. Is there any way toavoid storing all of theelements?Yes.
The user-installable output function capability is designedspecifically for this purpose. When you call the solver with nooutput arguments, the solver does not allocate storage to hold theentire solution history. Instead, the solver calls OutputFcn(t,y)at each time step. To keep the history of specific elements, writean output function that stores or plots only the elements you careabout.15-5115Differential EquationsProblem Size, Memory Use, and Computation Speed (Continued)QuestionAnswerWhat is the startup cost ofthe integration and howcan I reduce it?The biggest startup cost occurs as the solver attempts to find astep size appropriate to the scale of the problem. If you happen toknow an appropriate step size, use the InitialStep property. Forexample, if you repeatedly call the integrator in an event locationloop, the last step that was taken before the event is probably onscale for the next integration.
See ballode for an example.Time Steps for Integration15-52QuestionAnswerThe first step size that theintegrator takes is toolarge, and it missesimportant behavior.You can specify the first step size with the InitialStep property.The integrator tries this value, then reduces it if necessary.Can I integrate with fixedstep sizes?No.Initial Value Problems for ODEs and DAEsError Tolerance and Other OptionsQuestionAnswerHow do I choose RelTol andAbsTol?RelTol, the relative accuracy tolerance, controls the number ofcorrect digits in the answer.
AbsTol, the absolute error tolerance,controls the difference between the answer and the solution. Ateach step, the error e in component i of the solution satisfies|e(i)| <= max(RelTol*abs(y(i)),AbsTol(i))Roughly speaking, this means that you want RelTol correctdigits in all solution components except those smaller thanthresholds AbsTol(i). Even if you are not interested in acomponent y(i) when it is small, you may have to specifyAbsTol(i) small enough to get some correct digits in y(i) so thatyou can accurately compute more interesting components.I want answers that arecorrect to the precision ofthe computer.
Why can’t Isimply set RelTol to eps?You can get close to machine precision, but not that close. Thesolvers do not allow RelTol near eps because they try toapproximate a continuous function. At tolerances comparable toeps, the machine arithmetic causes all functions to lookdiscontinuous.How do I tell the solver thatI don’t care about gettingan accurate answer for oneof the solution components?You can increase the absolute error tolerance corresponding tothis solution component.
If the tolerance is bigger than thecomponent, this specifies no correct digits for the component. Thesolver may have to get some correct digits in this component tocompute other components accurately, but it generally handlesthis automatically.15-5315Differential EquationsSolving Different Kinds of ProblemsQuestionAnswerCan the solvers handlepartial differentialequations (PDEs) that havebeen discretized by themethod of lines?Yes, because the discretization produces a system of ODEs.Depending on the discretization, you might have a form involvingmass matrices – MATLAB ODE solvers provide for this.
Oftenthe system is stiff. This is to be expected when the PDE isparabolic and when there are phenomena that happen on verydifferent time scales such as a chemical reaction in a fluid flow. Insuch cases, use one of the four solvers: ode15s, ode23s, ode23t,ode23tb. If there are many equations, set the JPattern property.This might make the difference between success and failure dueto the computation being too expensive. When the system is notstiff, or not very stiff, ode23 or ode45 is more efficient thanode15s, ode23s, ode23t, or ode23tb.
For an example that usesJPattern, see “Example: Large, Stiff Sparse Problem” onpage 15-38.Parabolic-elliptic partial differential equations in 1-D can besolved directly with the MATLAB PDE solver, pdepe. See “PartialDifferential Equations” on page 15-82 for more information.15-54Can I solvedifferential-algebraicequation (DAE) systems?Yes. The solvers ode15s and ode23t can solve some DAEs of theform M ( t, y )y′ = f ( t, y ) where M ( t, y ) is singular.
The DAEsmust be of index 1. For examples, see amp1dae or hb1dae.Can I integrate a set ofsampled data?Not directly. You have to represent the data as a function byinterpolation or some other scheme for fitting data. Thesmoothness of this function is critical. A piecewise polynomial fitlike a spline can look smooth to the eye, but rough to a solver; thesolver takes small steps where the derivatives of the fit havejumps. Either use a smooth function to represent the data or useone of the lower order solvers (ode23, ode23s, ode23t, ode23tb)that is less sensitive to this.Initial Value Problems for ODEs and DAEsSolving Different Kinds of Problems (Continued)QuestionAnswerCan I solvedelay-differentialequations?Not directly. In some cases it is possible to use the initial valueproblem solvers to solve delay-differential equations by breakingthe simulation interval into smaller intervals the length of asingle delay.
For more information about this approach, see [6].What do I do when I havethe final and not the initialvalue?All the solvers of the MATLAB ODE suite allow you to solvebackwards or forwards in time. The syntax for the solvers is[t,y] = ode45(odefun,[t0 tf],y0);and the syntax accepts t0 > tf.TroubleshootingQuestionAnswerThe solution doesn’t looklike what I expected.If you’re right about its appearance, you need to reduce the errortolerances from their default values.
A smaller relative errortolerance is needed to compute accurately the solution ofproblems integrated over “long” intervals, as well as solutions ofproblems that are moderately unstable. You should checkwhether there are solution components that stay smaller thantheir absolute error tolerance for some time. If so, you are notasking for any correct digits in these components. This may beacceptable for these components, but failing to compute themaccurately may degrade the accuracy of other components thatdepend on them.My plots aren’t smoothenough.Increase the value of Refine from its default of 4 in ode45 and 1in the other solvers.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
