My IVP-Sheet 3 - CODING PDF

Preview

Summary

CODING...

Description

Solving an ordinary differential equation [≈25 lines of code] [See hints and example at the bottom of the description.] Write a Matlab function MyIVP that solves an initial-value problem (IVP) for a system of ordinary differential equations (ODEs) of the form x (t) = f (t, x (t)), ˙ where f : R × Rn 7→ Rn is an arbitrary function with one one-dimensional input (for time t ) and one n-dimensional input x, and n-dimensional output. The function should implement a Runge-Kutta formula (for example, the rk34 formula or the Dormand & Prince formula). The initial value x 0 is provided by the user of MyIVP. The first line of MyIVP (saved in a file MyIVP.m) should look like this function [xend,t,xt]=MyIVP(f,x0,tspan,N)

Inputs

• f: function defining the right-hand side of the ODE. f should accept two arguments: t (a number) and x (an n-dimensional vector). The function f should return an n-dimensional vector y (the time derivative). Typical calling sequence: y=f(t,x), returning the value of f at time t in position x. • x0: initial value where integration starts from (n-dimensional vector). • tspan: Starting time and end time for integration. Integration has to run from time t = tspan(1) to time t =tspan(2). • N: number of steps for integration. The integration stepsize h=(tspan(2)-tspan(1))/N should be small. Outputs

• xend: result of integration at t =tspan(2). • t: vector of times at which intermediate values have been computed (this should have N + 1 entries). • xt: intermediate values (n × (N + 1)-array). xt(:,k) should be the solution at t(k). You can check the built-in variable nargout inside your function to see if the user wants to get three outputs or only the end value xend. If nargout==1 you don’t need to store the intermediate values. Hints

First implement the skeleton of the function using the Euler formula x (t k+1 ) = x (t k ) + h f (t k , x (t k )) where t k =tspan(1)+k h. Test if the function works. Then try and improve on the Euler formula. Look up wikipedia for good Runge-Kutta formulas. The classical formula rk34 is decent, the Dormand & Prince pair is better. A better formula permits you to make larger steps, which will be helpful because some Coursework questions will require long computation times. Good formulas

Example

Suppose you want to solve the ODE ˙x 1 (t) = x 2 (t),

(1)

˙x 2 (t) = − x 1 (t) − x 2 (t) + cos t

with initial condition x ini = (x 1(0), x 2(0)) = (1, 0) for t from 0 to 2π . Here is a minimal script that shows how MyIVP can be used to solve this problem (note that “...” is used if one wants to continue a command on the next line) 1 2 3 4 5 6 7 8 9 10 11

%% Single test call of MyIVP clear %% x''=-x-x'+cos(t)

-> solution is [sin(t);cos(t)]

f=@(t,x)[x(2);...

% function f defining

-x(1)-x(2)+cos(t)];

% right-hand side

%% other inputs tspan=[0,2*pi]; N=20;

% number of steps

% time interval

xini=[0;1];

% initial value

[xend,t,xt]=MyIVP(f,xini,tspan,N); % call MyIVP plot(t,xt,'.-');

% plot solution

1

0.5

0

−0.5

−1 0

1

2

3

4

5

6

Figure 1: output of plot command in line 11 Try to call f (0, x ini ) in the command window after line 9 to see how the input f is supposed to work: >> f(0,xini) ans = 1 0 Testing All coursework questions explore differential equations, which makes a reliably working

function MyIVP important. For example system (1) we know that the solution is (x 1 (t ) , x 2(t)) = (sin t, cos t). Uploaded to ELE is a test script that runs your MyIVP for an increasing number of steps N over from time 0 to 2π and compares the difference to the known solution ( x end, yend = (0 , 1). A doubling of N divide the error by 32 for the Dormand & Prince Pair, by 16 for rk34, and by 2 for Euler....