Beginning Programming for Engineers

1. Beginning Programming for Engineers     Numerical Solutions to Differential Equations

2. 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

3. Euler's Method

4. 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.

5. 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);

6. 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).

7. 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);

8. 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 

9. 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.

10. 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.

11. 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);

12. 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].

13. 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)) ];

14. 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))

15. 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

16. 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.

17. 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

18. 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 n Event function to stop. • You can use more complicated function handles to pass information to the Event function, e.g., the terrain.

19. 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.