Using MATLAB (779505), страница 59
Текст из файла (страница 59)
The example assignsthis matrix to the property JPattern, and the solver uses the sparsity patternto generate the Jacobian numerically as a sparse matrix. Providing a sparsitypattern can significantly reduce the number of function evaluations requiredto generate the Jacobian and can accelerate integration. For the Brusselatorproblem, if the sparsity pattern is not supplied, 2 N evaluations of the functionare needed to compute the 2N-by-2N Jacobian matrix. If the sparsity pattern issupplied, only four evaluations are needed, regardless of the value of N.To run this example from the MATLAB Help browser, click on the examplename.
Otherwise, type brussode at the command line. At the command line,you can specify a value of N as an argument to brussode. The default isN = 20.function brussode(N)%BRUSSODE Stiff problem modeling a chemical reactionif nargin<1N = 20;endtspan = [0; 10];y0 = [1+sin((2*pi/(N+1))*(1:N)); repmat(3,1,N)];options = odeset('Vectorized','on','JPattern',jpattern(N));[t,y] = ode15s(@f,tspan,y0,options,N);u = y(:,1:2:end);x = (1:N)/(N+1);surf(x,t,u);view(-40,30);xlabel('space');ylabel('time');15-3915Differential Equationszlabel('solution u');title(['The Brusselator for N = ' num2str(N)]);% -------------------------------------------------------------function dydt = f(t,y,N)c = 0.02 * (N+1)^2;dydt = zeros(2*N,size(y,2));% preallocate dy/dt% Evaluate the two components of the function at one edge of% the grid (with edge conditions).i = 1;dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...c*(1-2*y(i,:)+y(i+2,:));dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...c*(3-2*y(i+1,:)+y(i+3,:));% Evaluate the two components of the function at all interior% grid points.i = 3:2:2*N-3;dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...c*(y(i-2,:)-2*y(i,:)+y(i+2,:));dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...c*(y(i-1,:)-2*y(i+1,:)+y(i+3,:));% Evaluate the two components of the function at the other edge% of the grid (with edge conditions).i = 2*N-1;dydt(i,:) = 1 + y(i+1,:).*y(i,:).^2 - 4*y(i,:) + ...c*(y(i-2,:)-2*y(i,:)+1);dydt(i+1,:) = 3*y(i,:) - y(i+1,:).*y(i,:).^2 + ...c*(y(i-1,:)-2*y(i+1,:)+3);% -------------------------------------------------------------function S = jpattern(N)B = ones(2*N,5);B(2:2:2*N,2) = zeros(N,1);B(1:2:2*N-1,4) = zeros(N,1);S = spdiags(B,-2:2,2*N,2*N);15-40Initial Value Problems for ODEs and DAEsThe Brusselator for N = 2032.5solution u21.510.5010180.860.640.42time0.200spaceExample: Simple Event Locationballode models the motion of a bouncing ball.
This example illustrates theevent location capabilities of the ODE solvers.The equations for the bouncing ball arey′1 = y2y′2 = – 9.8In this example, the event function is coded in a subfunction events[value,isterminal,direction] = events(t,y)which returns:• A value of the event function• The information whether or not the integration should stop (isterminal = 1or 0, respectively) when value = 0• The desired directionality of the zero crossings:15-4115Differential Equations-1Detect zero crossings in the negative direction only0Detect all zero crossings1Detect zero crossings in the positive direction onlyThe length of value, isterminal, and direction is the same as the number ofevent functions.
The ith element of each vector, corresponds to the ith eventfunction. For an example of more advanced event location, see orbitode(“Example: Advanced Event Location” on page 15-44).In ballode, setting the Events property to @events causes the solver to stop theintegration (isterminal = 1) when the ball hits the ground (the height y(1) is0) during its fall (direction = -1).
The example then restarts the integrationwith initial conditions corresponding to a ball that bounced.To run this example from the MATLAB Help browser, click on the examplename. Otherwise, type ballode at the command line.function ballode%BALLODE Run a demo of a bouncing ball.tstart = 0;tfinal = 30;y0 = [0; 20];refine = 4;options = odeset('Events',@events,'OutputFcn', @odeplot,...'OutputSel',1,'Refine',refine);set(gca,'xlim',[0 30],'ylim',[0 25]);box onhold on;tout = tstart;yout = y0.';teout = [];yeout = [];ieout = [];for i = 1:10% Solve until the first terminal event.[t,y,te,ye,ie] = ode23(@f,[tstart tfinal],y0,options);15-42Initial Value Problems for ODEs and DAEsif ~isholdhold onend% Accumulate output.nt = length(t);tout = [tout; t(2:nt)];yout = [yout; y(2:nt,:)];teout = [teout; te];% Events at tstart are never reported.yeout = [yeout; ye];ieout = [ieout; ie];ud = get(gcf,'UserData');if ud.stopbreak;end% Set the new initial conditions, with .9 attenuation.y0(1) = 0;y0(2) = -.9*y(nt,2);% A good guess of a valid first time step is the length of% the last valid time step, so use it for faster computation.options = odeset(options,'InitialStep',t(nt)-t(nt-refine),...'MaxStep',t(nt)-t(1));tstart = t(nt);endplot(teout,yeout(:,1),'ro')xlabel('time');ylabel('height');title('Ball trajectory and the events');hold offodeplot([],[],'done');% -------------------------------------------------------------function dydt = f(t,y)dydt = [y(2); -9.8];% -------------------------------------------------------------function [value,isterminal,direction] = events(t,y)% Locate the time when height passes through zero in a% decreasing direction and stop integration.15-4315Differential Equationsvalue = y(1);isterminal = 1;direction = -1;% Detect height = 0% Stop the integration% Negative direction onlyBall trajectory and the events2520height151050051015time202530Example: Advanced Event Locationorbitode illustrates the solution of a standard test problem for those solversthat are intended for nonstiff problems.
It traces the path of a spaceshiptraveling around the moon and returning to the earth. (Shampine andGordon [7], p.246).The orbitode problem is a system of the four equations shown below15-44Initial Value Problems for ODEs and DAEsy′1 = y 3y′2 = y 4µ∗ ( y 1 + µ ) µ ( y 1 – µ∗ )y′3 = 2y 4 + y 1 – --------------------------- – --------------------------r 23r 31µ∗ y 2 µy 2y′4 = – 2y 3 + y 2 – ------------ – --------r 31r 23whereµ = 1 ⁄ 82.45µ∗ = 1 – µr1 =( y 1 + µ ) 2 + y 22r2 =( y 1 – µ∗ ) + y 222The first two solution components are coordinates of the body of infinitesimalmass, so plotting one against the other gives the orbit of the body. The initialconditions have been chosen to make the orbit periodic.
The value of µcorresponds to a spaceship traveling around the moon and the earth.Moderately stringent tolerances are necessary to reproduce the qualitativebehavior of the orbit. Suitable values are 1e-5 for RelTol and 1e-4 for AbsTol.The events subfunction includes event functions which locate the point ofmaximum distance from the starting point and the time the spaceship returnsto the starting point. Note that the events are located accurately, even thoughthe step sizes used by the integrator are not determined by the location of theevents. In this example, the ability to specify the direction of the zero crossingis critical. Both the point of return to the initial point and the point ofmaximum distance have the same event function value, and the direction of thecrossing is used to distinguish them.To run this example from the MATLAB Help browser, click on the examplename.
Otherwise, type orbitode at the command line. The example uses the15-4515Differential Equationsoutput function odephase2 to produce the two-dimensional phase plane plotand let you to see the progress of the integration.function orbitode%ORBITODE Restricted three-body problemtspan = [0 7];y0 = [1.2; 0; 0; -1.04935750983031990726];options = odeset('RelTol',1e-5,'AbsTol',1e-4,...'OutputFcn',@odephas2,'Events',@events);[t,y,te,ye,ie] = ode45(@f,tspan,y0,options);plot(y(:,1),y(:,2),ye(:,1),ye(:,2),'o');title ('Restricted three body problem')ylabel ('y(t)')xlabel ('x(t)')% -------------------------------------------------------------function dydt = f(t,y)mu = 1 / 82.45;mustar = 1 - mu;r13 = ((y(1) + mu)^2 + y(2)^2) ^ 1.5;r23 = ((y(1) - mustar)^2 + y(2)^2) ^ 1.5;dydt = [ y(3)y(4)2*y(4) + y(1) - mustar*((y(1)+mu)/r13) - ...mu*((y(1)-mustar)/r23)-2*y(3) + y(2) - mustar*(y(2)/r13) - mu*(y(2)/r23) ];% -------------------------------------------------------------function [value,isterminal,direction] = events(t,y)% Locate the time when the object returns closest to the% initial point y0 and starts to move away, and stop integration.% Also locate the time when the object is farthest from the% initial point y0 and starts to move closer.%% The current distance of the body is%%DSQ = (y(1)-y0(1))^2 + (y(2)-y0(2))^2%= <y(1:2)-y0,y(1:2)-y0>%15-46Initial Value Problems for ODEs and DAEs% A local minimum of DSQ occurs when d/dt DSQ crosses zero% heading in the positive direction.
We can compute d(DSQ)/dt as%% d(DSQ)/dt = 2*(y(1:2)-y0)'*dy(1:2)/dt = 2*(y(1:2)-y0)'*y(3:4)%y0 = [1.2; 0];dDSQdt = 2 * ((y(1:2)-y0)' * y(3:4));value = [dDSQdt; dDSQdt];isterminal = [1; 0];% Stop at local minimumdirection = [1; -1];% [local minimum, local maximum]Restricted three body problem0.80.60.4y(t)0.20−0.2−0.4−0.6−0.8−1.5−1−0.50x(t)0.511.5Example: Differential-Algebraic Problemhb1dae reformulates the hb1ode example as a differential-algebraic equation(DAE) problem.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
