Using MATLAB (779505), страница 58
Текст из файла (страница 58)
The default value of maybe causes thesolver to test whether the problem is a DAE, i.e.,whether M(t 0, y 0) is singular.For an example of a problem with a mass matrix,see “Example: Finite Element Discretization” onpage 15-35 (fem1ode).InitialSlopeVector |{zero vector}Vector representing the consistent initial slope yp 0 ,where yp 0 satisfies M(t 0, y 0) yp 0 = f(t 0, y 0) . Thedefault is the zero vector.This property is for use with the ode15s and ode23tsolvers when solving DAEs.Event Location PropertyIn some ODE problems the times of specific events are important, such as thetime at which a ball hits the ground, or the time at which a spaceship returnsto the earth.
While solving a problem, the MATLAB ODE solvers can detectsuch events by locating transitions to, from, or through zeros of user-definedfunctions.15-2915Differential EquationsThe following table describes the Events property. Use odeset to set thisproperty.ODE Events PropertyStringValueDescriptionEventsFunctionMATLAB function that includes one or more event functions. TheMATLAB function is of the form[value,isterminal,direction] = events(t,y)value, isterminal, and direction are vectors for which the ithelement corresponds to the ith event function:• value(i) is the value of the ith event function.• isterminal(i) = 1 if the integration is to terminate at a zero ofthis event function and 0 otherwise.• direction(i) = 0 if all zeros are to be located (the default), +1 ifonly zeros where the event function is increasing, and -1 if onlyzeros where the event function is decreasing.When you specify an events function, the solver returns threeadditional outputs, TE, YE, and IE.options = odeset(‘Events’,@events)[T,Y,TE,YE,IE] = solver(odefun,tspan,y0,options)• TE is a column vector of times at which events occur.• Rows of YE are the solution values corresponding to times in TE.• Indices in vector IE specify which event occurred at the time in TE.For examples that use an event function, see “Example: SimpleEvent Location” on page 15-41 (ballode) and “Example: AdvancedEvent Location” on page 15-44 (orbitode).ode15s Propertiesode15s is a variable-order solver for stiff problems.
It is based on the numericaldifferentiation formulas (NDFs). The NDFs are generally more efficient thanthe closely related family of backward differentiation formulas (BDFs), also15-30Initial Value Problems for ODEs and DAEsknown as Gear’s methods. The ode15s properties let you choose between theseformulas, as well as specifying the maximum order for the formula used.The following table describes the ode15s properties. Use odeset to set theseproperties.ode15s PropertiesPropertyValueDescriptionMaxOrder1 | 2 | 3 | 4 | {5}The maximum order formula used to compute the solution.BDFon | {off}Specifies whether you want to use the BDFs instead of thedefault NDFs.
Set BDF on to have ode15s use the BDFs.For both the NDFs and BDFs, the formulas of orders 1 and2 are A-stable (the stability region includes the entire lefthalf complex plane). The higher order formulas are not asstable, and the higher the order the worse the stability.There is a class of stiff problems (stiff oscillatory) that issolved more efficiently if MaxOrder is reduced (for exampleto 2) so that only the most stable formulas are used.Examples: Applying the ODE Initial Value ProblemSolversThis section contains several examples that illustrate the kinds of problemsyou can solve in MATLAB:• Simple nonstiff problem (rigidode)• Stiff problem (vdpode)• Finite element discretization (fem1ode)• Large, stiff sparse problem (brussode)• Simple event location (ballode)• Advanced event location (orbitode)• Differential-algebraic problem (hb1dae)15-3115Differential EquationsExample: Simple Nonstiff Problemrigidode illustrates the solution of a standard test problem proposed by Kroghfor solvers intended for nonstiff problems [7].The ODEs are the Euler equations of a rigid body without external forces.y′1 = y2 y 3y′2 = – y 1 y 3y′3 = – 0.51 y 1 y 2For your convenience, the entire problem is defined and solved in a singleM-file.
The differential equations are coded as a subfunction f. Because theexample calls the ode45 solver without output arguments, the solver uses thedefault output function odeplot to plot the solution components.To run this example from the MATLAB Help browser, click on the examplename. Otherwise, type rigidode at the command line.function rigidode%RIGIDODE Euler equations of a rigid body without external forcestspan = [0 12];y0 = [0; 1; 1];% Solve the problem using ode45ode45(@f,tspan,y0);% -----------------------------------------------------------function dydt = f(t,y)dydt = [ y(2)*y(3)-y(1)*y(3)-0.51*y(1)*y(2) ];15-32Initial Value Problems for ODEs and DAEs10.80.60.40.20−0.2−0.4−0.6−0.8−1024681012Example: Stiff Problem (van der Pol Equation)vdpode illustrates the solution of the van der Pol problem described in“Example: The van der Pol Equation, µ = 1000 (Stiff)” on page 15-14.
Thedifferential equationsy′1 = y22y′2 = µ(1 – y 1 )y 2 – y 1involve a constant parameter µ .As µ increases, the problem becomes more stiff, and the period of oscillationbecomes larger. When µ is 1000 the equation is in relaxation oscillation andthe problem is very stiff.
The limit cycle has portions where the solutioncomponents change slowly and the problem is quite stiff, alternating withregions of very sharp change where it is not stiff (quasi-discontinuities).By default, the solvers in the ODE suite that are intended for stiff problemsapproximate Jacobian matrices numerically. However, this example providesa subfunction J(t,y,mu) to evaluate the Jacobian matrix ∂f ⁄ ∂y analytically at15-3315Differential Equations(t,y) for µ = mu. The use of an analytic Jacobian can improve the reliabilityand efficiency of integration.To run this example from the MATLAB Help browser, click on the examplename.
Otherwise, type vdpode at the command line. At the command line, youcan specify a value of µ as an argument to vdpode. The default is µ = 1000.function vdpode(MU)%VDPODE Parameterizable van der Pol equation (stiff for large MU)if nargin < 1MU = 1000;% defaultendtspan = [0; max(20,3*MU)];y0 = [2; 0];options = odeset('Jacobian',@J);% Several periods[t,y] = ode15s(@f,tspan,y0,options,MU);plot(t,y(:,1));title(['Solution of van der Pol Equation, \mu = ' num2str(MU)]);xlabel('time t');ylabel('solution y_1');axis([tspan(1) tspan(end) -2.5 2.5]);--------------------------------------------------------------function dydt = f(t,y,mu)dydt = [y(2)mu*(1-y(1)^2)*y(2)-y(1) ];--------------------------------------------------------------function dfdy = J(t,y,mu)dfdy = [01-2*mu*y(1)*y(2)-1mu*(1-y(1)^2) ];15-34Initial Value Problems for ODEs and DAEsSolution of van der Pol Equation, µ = 10002.521.51solution y10.50−0.5−1−1.5−2−2.5050010001500time t200025003000Example: Finite Element Discretizationfem1ode illustrates the solution of ODEs that result from a finite elementdiscretization of a partial differential equation.
The value of N in the callfem1ode(N) controls the discretization, and the resulting system consists of Nequations. By default, N is 19.This example involves a mass matrix. The system of ODEs comes from amethod of lines solution of the partial differential equation2– t ∂u∂ue ------ = --------∂t ∂x 2with initial condition u ( 0, x ) = sin ( x ) and boundary conditionsu ( t, 0 ) = u ( t, π ) = 0 . An integer N is chosen, h is defined as π ⁄ ( N + 1 ) , andthe solution of the partial differential equation is approximated at x k = kh fork = 0, 1, ..., N+1 by15-3515Differential EquationsNu (t,x k) ≈∑ ck ( t )φk ( x )k=1Here φ k(x) is a piecewise linear function that is 1 at x k and 0 at all the otherx j .
A Galerkin discretization leads to the system of ODEsM ( t )c′ = Jc where c ( t ) =c1 ( t )cN ( t )and the tridiagonal matrices M ( t ) and J are given by 2h ⁄ 3 exp ( – t )M ij = h ⁄ 6 exp ( – t )0if i = jif i = j ± 1otherwiseand –2 ⁄ hJ ij = 1 ⁄ h 0if i = jif i = j ± 1otherwiseThe initial values c ( 0 ) are taken from the initial condition for the partialdifferential equation. The problem is solved on the time interval [ 0, π ] .In fem1ode the propertiesoptions = odeset('Mass',@mass,'MStateDep','none','Jacobian',J)indicate that the problem is of the form M(t)y′ = Jy . The subfunctionmass(t,N) evaluates the time-dependent mass matrix M ( t ) and J is theconstant Jacobian.To run this example from the MATLAB Help browser, click on the examplename. Otherwise, type fem1ode at the command line.
At the command line, youcan specify a value of N as an argument to fem1ode. The default is N = 19.function fem1ode(N)%FEM1ODE Stiff problem with a time-dependent mass matrix15-36Initial Value Problems for ODEs and DAEsif nargin < 1N = 19;endh = pi/(N+1);y0 = sin(h*(1:N)');tspan = [0; pi];%edJThe Jacobian is constant.= repmat(1/h,N,1);% e=[(1/h) ... (1/h)];= repmat(-2/h,N,1);% d=[(-2/h) ... (-2/h)];= spdiags([e d e], -1:1, N, N);options = odeset('Mass',@mass,'MStateDependence','none', ...'Jacobian',J);[t,y] = ode15s(@f,tspan,y0,options,N);surf((1:N)/(N+1),t,y);set(gca,'ZLim',[0 1]);view(142.5,30);title(['Finite element problem with time-dependent mass ' ...'matrix, solved by ODE15S']);xlabel('space ( x/\pi )');ylabel('time');zlabel('solution');%--------------------------------------------------------------function out = f(t,y,N)h = pi/(N+1);e = repmat(1/h,N,1);% e=[(1/h) ...
(1/h)];d = repmat(-2/h,N,1);% d=[(-2/h) ... (-2/h)];J = spdiags([e d e], -1:1, N, N);out = J*y;%--------------------------------------------------------------function M = mass(t,N)h = pi/(N+1);e = repmat(exp(-t)*h/6,N,1); % e(i)=exp(-t)*h/6e4 = repmat(4*exp(-t)*h/6,N,1);M = spdiags([e e4 e], -1:1, N, N);15-3715Differential EquationsFinite element problem with time−dependent mass matrix, solved by ODE15S10.8solution0.60.40.200010.220.40.630.841space ( x/π )timeExample: Large, Stiff Sparse Problembrussode illustrates the solution of a (potentially) large stiff sparse problem.The problem is the classic “Brusselator” system [3] that models diffusion in achemical reaction22u′i = 1 + u i v i – 4u i + α ( N + 1 ) ( u i – 1 – 2u i + u i + 1 )22v′i = 3u i – u i v i + α ( N + 1 ) ( v i – 1 – 2v i + v i + 1 )and is solved on the time interval [0,10] with α = 1/50 andu i ( 0 ) = 1 + sin ( 2πx i ) vi ( 0 ) = 3with x i = i ⁄ ( N + 1 ), for i = 1, ..., NThere are 2N equations in the system, but the Jacobian is banded with aconstant width 5 if the equations are ordered as u 1, v 1, u 2, v 2, …15-38Initial Value Problems for ODEs and DAEsIn the call brussode(N), where N corresponds to N , the parameter N ≥ 2specifies the number of grid points.
The resulting system consists of 2Nequations. By default, N is 20. The problem becomes increasingly stiff and theJacobian increasingly sparse as N increases.The subfunction f(t,y,N) returns the derivatives vector for the Brusselatorproblem. The subfunction jpattern(N) returns a sparse matrix of 1s and 0sshowing the locations of nonzeros in the Jacobian ∂f ⁄ ∂y .
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
