1 / 20

# Beginning Programming for Engineers

Numerical Solutions to Differential Equations. Beginning Programming for Engineers    . Learning goals for class 5. Understand why numeric solutions of ODEs is necessary.

Télécharger la présentation

## Beginning Programming for Engineers

An Image/Link below is provided (as is) to download presentation Content is provided to you AS IS for your information and personal use only. Download presentation by click this link. During download, if you can't get a presentation, the file might be deleted by the publisher.

E N D

### Presentation Transcript

1. Numerical Solutions to Differential Equations Beginning Programming for Engineers

2. Learning goals for class 5 • Understand why numeric solutions of ODEs is necessary. • Understand, by referring to Euler's method, how it could be possible to numerically solve an ODE, even when their is no algebraic solution. • Know how to use Matlab's ODE solvers to solve ODE IVP problems, with one or more equations, and with higher level derivatives; in particular, how to use ode45. • Understand why more than one ODE solver is needed, and have some advice how to choose.

3. Initial Value Problems We wish to "solve" the ordinary differential equation (ODE): where we have the initial conditions: By "solve", we want to compute            for

4. Euler's Method

5. Matlab's ODE Solvers Matlab provides a selection of solvers to integrate an ODE function over a range.  General form: [t,y] = solver(odefun, tspan, y0) The odefun is the function to be integrated over the interval tspan with initial conditions y0.  The result is the vector t for values of the independent variable in the range tspan, and estimated values y for each value t. The solver is typically ode45.

6. Simple usage of ode45 We would like to solve this system numerically: We need a function to calculate function [dy] = simple_ode(t, y)    dy = -3.7*y;end Call the solver: [t,y] = ode45(@simple_ode, [0 1], 5);

7. Using ode45, continued Given the example: [t,y] = ode45(@simple_ode, [0 1], 5); • The interval [0 1] is the range over which the independent variable (often x or t) varies.  The solver chooses its own values to use.  Actual values chosen are sorted in t, a column vector. • Corresponding values of the integrated function are stored in the column vector y. • Instead of [0 1], additional points can be given, e.g, 0:0.1:1 -- which will be used as the values of the output t (and to compute y).

8. Function Handles The notation @simple_ode is a function handle, to allow the solver to call our function simple_ode. Small functions can be defined anonymously: odefun = @(T,Y) -3.7*Y; [t,y] = ode45(odefun, [0 1], 5); or even: [t,y] = ode45(@(T,Y)-3.7*Y, [0 1], 5);

9. Multiple ODEs The Matlab ODE solvers can solve systems of related equations over the same interval of the independent variable. Example: The Lotka-Volterra predator-prey model can be used to model the population on rabbits (r) and foxes (f): Model this with

10. Rabbits and Foxes (2) ... Need the ODE function: function [dy] = foxrabbitode(t,y)     % FOXRABBITODE calculates population changes     % y(1) = rabbits     % y(2) = foxes     alpha = 0.01;     r = y(1);     f = y(2);     dy = [ 2*r - alpha*r*f            -f + alpha*r*f ]; end • Calculate on y, a column vector, and produce a column vector.  Each row corresponds to one of the equations.

11. Rabbits and Foxes (3) ... Create column vector of initial conditions (r0 and f0) y0 = [ 300; 150 ]; [t,y] = ode45(@foxrabbitode, [0 20], y0); figure; hold on; plot(t, y(:,1), 'b');  % rabbits plot(t, y(:,2), 'r');  % foxes Each column of dy corresponds to the points computed for one of the equations.

12. Function handles and arguments We might want to pass alpha as an argument to the ODE function: function [dy] = foxrabbitode2(t,y,alpha)     r = y(1);     f = y(2);     dy = [ 2*r - alpha*r*f            -f + alpha*r*f ]; end Need anonymous function to pass this to ODE function: odefun = @(T,Y) foxrabbitode2(T,Y,alpha);[t,y] = ode45(odefun, interval, y0);

13. Higher-order ODEs Treat higher-order ODEs as systems of related ODEs. The ODE function receives both y and y', and computes y' and y''.  (Simply copy the received y' value!) We need initial conditions for y and y'. Consider solving: Let: y(1) = 0,  y'(1) = 1 Solve for y(t), for t in [1 30].

14. Function for simple second-order ODE function [dy] = second_ode(t,y)     % SECOND_ODE is a simple second order ODE.     %     % y(1) = f(t)    % 1st equation     % y(2) = f'(t)   % 2nd equation     %     % dy(1) = f'(t)  % 1st equation     % dy(2) = f''(t) % 2nd equation     % Note how we can just copy y(2) to dy(1).     dy = [ y(2)           (-1/t)*sin(y(1)) ];

15. Solving the second-order ODE % Compute the solutions. [t,y] = ode45(@second_ode, [1 30], [0;1]); % y(:,1) is y(t), and y(:,2) is y'(t). plot(t, y(:,1))

16. Trajectory as ODE function [dpv] = update_path(t, pv)   %    pv = position and velocity state vector  %        pv(1) = x position, xpos, at time t  %        pv(2) = y position, ypos, at time t  %        pv(3) = x velocity, vx, at time t  %        pv(4) = y velocity, vy, at time t   x = pv(1);   y = pv(2);   xv = pv(3);   yv = pv(4);   dpv = [xv        % dx = x velocity         yv        % dy = y velocity         0         % d2x = d(x velocity) = x acceleration         -9.81 ];  % d2y = d(y velocity) = y acceleration

17. Event functions Sometimes we need to stop the integration due to some event, especially if we don't know how long to run! You can create an "event function" to stop integration: function [value,isterminal,direction] = event(t,y) This accepts t,y like the ODE function.  It computes a value at t,y.  The output isterminal, if 1, means to stop the integration if value is 0.  The direction indicates when we care about value being 0; if it is 0, then we always care.

18. Event functions: Trajectory Sample for testing if above ground: function [value, isterminal, direction] = ...          trajectory_event(t, pv)       if pv(2) > 0        value = 1;  % Above the ground    else        value = 0;  % Hit the ground    end    isterminal = 1; % Stop when value=0    direction = 0;  % Always care about zerosend

19. Using an event function Use odeset to configure options to the ODE solvers, including events functions: options = odeset('Event', @trajectory_event);[t,y] = ode45(@update_path, 0:dt:ft, ...                   [x;y;v0x;v0y], options); • You can use [0 Inf] as an interval when you don't know when to finish, but have an Event function to stop. • You can use more complicated function handles to pass information to the Event function, e.g., the terrain.

20. Other solvers • Some problems are stiff, needing slower specialized solvers.  Other solvers are faster than ode45, but not as accurate. • Usually try ode45 first.  If it works but is too slow, consider ode23. • If ode45 is so slow it is unusable, you might have a stiff problem.  Try using ode15s or ode23s.  • See doc ode45 for much more information, including guidance about choosing a solver, options you can pass to solvers, etc.

More Related