LAB 2

Table of Contents

2. Begynnelsevärdesproblem 2.1 Riktningsfält 2.2 Differensmetoder 2.3 Konvergens Uppgift 2 Uppgift 3: Jämför euler framåtmetod med ode45(). 2.5 System av ODE med ode45() Fasporträtt: u_1 vs. u_2 Limit cycle Uppgift 4: same procedure for new system of ODE's. 2.6: Högre ordningens diferentialekvationer Uppgift 6: Damped pendulum Randvärdesproblem/Boundary-value problems 3.2 Inskjutningsmetoden/Shooting method (https://en.wikipedia.org/wiki/Shooting_method) Solve the previous heat transfer problem using the Shooting method

2. Begynnelsevärdesproblem

Initial-value problems for first order differential equations have the form

, where u0 is a constant.

What makes it a differential equation is that u' is a function of u. The function will evolve differently depending on the initial value.

2.1 Riktningsfält

f=@(t,u) -u+sin(t)+cos(t); % högerled. dvs. u'(t)
f_exakt=@(t,u0) sin(t)+u0.*exp(-t); % Analytic uprim.

Draw vector-field of

. This shows how the function will change at each point.

That is, the axes are u(t) and t.

close all

a=0;

b=4;

N=20;

t=linspace(a,b,N);

u=linspace(-2,2,N); % u(t)

[T,U]=meshgrid(t,u); % Grid of (t, u(t))

DT=ones(size(T)); % d/dt(t) = 1

DU=f(T,U); % For our range of t and u(t), find u'(t).

figure

quiver(T,U,DT,DU,0.9) % Plot vector [dt du] at each point.

hold on

xlabel('t');

ylabel('u(t)')

hold on

grid on

2.2 Differensmetoder

We approximate the derivative at a point using a forward difference quotient: a linearisation between two adjacent elements.

, an approximation of the definition given at the top of the script.

Rearranging gives

If we let

be an approximation of

, then Euler's forward-method is

This method is explicit since all variables on the r.h.s. are known.

We derive using Euler's forward-method, using initial value

, and plot the curve.

N = 20;
h = (b-a)/N;
t = linspace(a,b,N+1);
u = zeros(size(t));
u0 = 1;
u(1) = u0; % Initial value.
for n = 1:N
    u(n+1) = u(n)+h.*f(t(n), u(n));
end
plot(t, u);
xlabel('t');
ylabel('u(t)')
grid on

Euler's backwards-method is derived in a similar way:

This method is implicit because

is unknown and appears on both sides!

In each step, a non-linear equation

must be solved using eg. newton's method.

The trapezium method combines Euler's backwards and forwards methods, and is thus also implicit.

Heun's Method takes the trapezium method and substitutes the implicit

for a forward method approximation (and substitutes

), making it explicit.

u = zeros(size(t));

u(1) = u0;

for n = 1:N

u(n+1) = u(n) + h/2*( f(t(n),u(n)) + f(t(n)+h, u(n)+h*f(t(n),u(n))) );

end

plot(t, u)

legend("","Euler forward", "Heun")

title(['u_0 = ' num2str(u0)])

%plot(t, f(t, 0), 'g--',t, f(t, -1), 'g--',t, f(t, 1), 'r-',t, f(t, 2), 'g--' )

%hold off

hold off

2.3 Konvergens

All the difference methods described above converge, and thus we can get arbitrarily close to the exact value as

Euler's methods have linear convergence, while Heun's method has quadratic convergence.

That is, if we halve the step-size, the error will be 1/2 and 1/4, respectively.

Uppgift 2

Lös problemet med Eulers framåtekvation. Välj en bra steglängd. Rita riktningsfältet.

close all
uprim=@(t,u) cos(3.*t)-sin(5*t).*u;
a=0;
b=5;
N=20;
t=linspace(a,b,N); % För N intervall krävs N+1 punkter.
u=linspace(0, 3, N);
[T,U]=meshgrid(t,u);
DT=ones(size(T)); DU=uprim(T,U);
quiver(T,U,DT,DU, 0.8), 
hold on

Vi ser att Euler metoden konvergerar mot den exakta lösningen i takt med att steglängden minskas.

Vi kan även se att felet ackumuleras med tid (dvs. med varje nytt steg i diff.ekvationen).

for N=[10 20 30 50 100 1000 10000]

h=(b-a)/N;

t=linspace(a,b,N+1);u=zeros(size(t));

ua=2;

u(1)=ua; % at u(t), t=0

for n=1:N

u(n+1)=u(n)+h*uprim(t(n), u(n));

end

plot(t, u)

end

grid on

xlabel('t');

ylabel('u(t)')

hold off

We can also look at the resulting value of the first equation

Depending on the method and step-size used, different errors will be introduced.

We will look at the difference between forward-Euler, Heun and the analytic solution at t=4.

close all

u0 = 1;

u_exakt = @(t,u0) sin(t) + u0*exp(-t);

uprim = @(t,u) -u +sin(t) + cos(t);

a = 0;

b = 4;

h_vec = [0.4 0.2 0.1 0.05 0.025 0.01 0.001];

N = (b-a)./h_vec;

res_eulerForward = zeros(size(N));

res_heun = zeros(size(N));

res_exakt = 0

res_exakt = 0

for i = 1:length(N)

t = linspace(a,b,N(i)+1);

u = zeros(size(t));

h = (b-a)./N(i);

u(1) = u0;

% Euler forward

for n=1:N(i)

u(n+1)=u(n)+h*uprim(t(n), u(n));

end

%plot(t, u); hold on

res_eulerForward(i) = u(end);

u = zeros(size(t));

u(1) = u0;

for n = 1:N(i)

u(n+1) = u(n) + h/2*( uprim(t(n),u(n)) + uprim(t(n)+h, u(n)+h*uprim(t(n),u(n))) );

end

res_heun(i) = u(end);

end

error_eulerForward = res_eulerForward - u_exakt(4,u0);

plot(h_vec, error_eulerForward); hold on

error_heun = res_heun - u_exakt(4,u0);

plot(h_vec, error_heun)

grid on

title("Convergence of forward Euler and Heun")

xlabel("h")

set ( gca, 'xdir', 'reverse' ) % Flip x axis in plot.

legend('Euler Forward', 'Heun')

Uppgift 3: Jämför euler framåtmetod med ode45().

close all

figure

uprim=@(t,u) cos(3.*t)-sin(5*t).*u;

a=0;

b=5;

N=20;

t=linspace(a,b,N+1); % För N intervall krävs N+1 punkter.

u=linspace(0, 3, N);

[T,U]=meshgrid(t,u);

DT=ones(size(T)); DU=uprim(T,U);

quiver(T,U,DT,DU, 0.8),

hold on

h=(b-a)/N;

u=zeros(size(t));

u(1)=ua; % at u(t), t=0

for n=1:N

u(n+1)=u(n)+h*uprim(t(n), u(n));

end

plot(t, u), hold on

[t,u]=ode45(uprim,[a b], ua);

plot(t,u), grid on

legend('','Eulers framåtmetod', 'ode45')

2.5 System av ODE med ode45()

Previously, our definition of an ODE was

This can be extended to systems of ODE's without changing the definition, by using vector notation.

The above says:

is a system of equations

that we are trying to obtain/numerically compute.

is a system of differential equations that are functions of t and the equations in

are the initial values of the equations defined in

Ex.

We have the following system of diff. equations:

After using ode45(), we obtain numerical approximations of

, the original functions.

close all

f1 = @(u1,u2) 2*u1-u1.*u2; % The differential equations.

f2 = @(u1,u2) u2+0.4*u1.*u2-u2.^2;

f = @(t,u) [f1(u(1),u(2)); f2(u(1),u(2))]; % Vectorized. OBS! odefun must accept arguments (t,y), even if one of them is not used!

u0 = [0.1; 0.1]; % initial values.

T = 8; % Max time bound.

tspan = [0 T];

[t, U] = ode45(f, tspan, u0);

plot(t,U(:,1), 'b', t, U(:,2), 'g'), hold on, grid on

hold off

legend('u_1(t)', 'u_2(t)')

This system of ODEs is autonomous, that is

, since the derivative does not directly depend on the indpendent variable t, only indirectly through u(t). That means that the evolution of the system depends only on the current state of the system.

This means that stationary points can be found. The stationary points are shown in the phase portrait.

If the independent variable is time, the system is said to be time-invariant.

Fasporträtt: u_1 vs. u_2

The phase portrait can be interpreted as visualizing the trajectories of solutions to the system depending on initial conditions for

In the phase portrait below, we see three stationary points. One is stable, corresponding to the final, stable solution of the time-plot above. Some initial conditions will evolve to this stable point, others will diverge.

u1 = linspace(-1,4,25);

u2 = linspace(-1,3,25);

[U1,U2] = meshgrid(u1,u2);

quiver(U1,U2,f1(U1,U2),f2(U1,U2), 1.5), hold on, grid on % vector components are the derivatives u1' and u2'.

plot(U(:,1), U(:,2), 'r') % Plot the solution that corresponds to a particular initial condition u0.

% To illustrate, change initial condition and recompute the system of ODEs.

% A different trajectory will be followed.

% It stabilizes at the same point (sink).

u0 = [0.1 0.5 -0.5 2; 0.5 3 3 3]; % Iterate over these pairs of points (columns are the pairs).

for i = 1:length(u0)

[t, U] = ode45(f, tspan, u0(:,i));

plot(U(:,1), U(:,2), 'g-', 'LineWidth',2)

end

xlabel('u_1'); ylabel('u_2');

axis([-1 4 -1 3])

plot(0,0,'r.', 0,1,'r.', 2.5 ,2,'r.')

legend('', 'u_0=[0.1;0.1]', 'u_0=[0.1;0.5]')

Limit cycle

As another example, we solve another autonomous that results in a single unstable stationary point.

Solutions (based on initial conditions) will follow a trajectory that spirals into a limit cycle (a closed trajectory that other solutions spiral into as

The stationary point resides in the interior of the limit cycle.

close all

f1 = @(u1, u2) -u1 + 0.04*u2 + u1.^2.*u2;

f2 = @(u1, u2) 0.8 - 0.04*u2 - u1.^2.*u2;

f = @(t,u) [f1(u(1),u(2)); f2(u(1),u(2))]; % System of diff.equations.

u1 = linspace(0,4.5,20);

u2 = linspace(0, 3, 20);

[U1, U2] = meshgrid(u1, u2);

quiver(U1, U2, f1(U1, U2), f2(U1, U2),5.5) % How to normalize arrow length? ie. scale each arrow

hold on

tspan = [0 12];

u0 = [0.1 0.5 1 2 1 ; 0.5 3 3 3 1.5]; % Iterate over these pairs of points (columns are the pairs).

for i = 1:length(u0)

[t, U] = ode45(f, tspan, u0(:,i));

plot(u0(1,i), u0(2,i), 'ko')

plot(U(:,1), U(:,2), 'LineWidth',2)

end

grid on

plot(0.8,1.1765, 'r.', 'LineWidth', 2)

axis([0 4.5 0 3])

Uppgift 4: same procedure for new system of ODE's.

This system has two stable stationary points.

% Funktioner u'1 och u'2 i ekvationssystem

f1=@(u1,u2) cos(1+u1-u2.^2+0.3);

f2=@(u1,u2) sin(0.3+u1+u2-0.4*u1);

f=@(t,u) [f1(u(1),u(2)); f2(u(1),u(2))];

% Begynnelsevärden u1(0) och u2(0)

u0=[1.5; 3];

% Lösningsintervall

T=8; tspan=linspace(0,T,200);

% Lös ekv.systemet med ode45 på interval tspan med begynnelsevärden u0.

[t,U]=ode45(f,tspan,u0);

% Tidsfunktion

subplot(1,3,1)

plot(t,U(:,1),t,U(:,2)), grid on, hold on, legend('u_1', 'u_2'), xlabel('t')

% Riktiningsfält

% First, set up solution space

u1=linspace(0,3,30);

u2=linspace(0,3,30);

[U1,U2]=meshgrid(u1,u2);

% Then use quiver with the individual functions, and our grid as points, as

% parameters.

subplot(1,3,2:3)

quiver(U1,U2,f1(U1,U2),f2(U1,U2), 0.9), hold on

axis([0 3 0 3]), axis square, xlabel('u_1'), ylabel('u_2');

% Riktningsbanor

for n1 = 0:1:3

for n2 = 0:1:3

u0(1)=n1;

u0(2)=n2;

[t,U]=ode45(f,tspan,u0);

plot(U(:,1),U(:,2), 'r', 'LineWidth', 1.5)

end

grid on

2.6: Högre ordningens diferentialekvationer

We can rewrite a higher order diff.equation as a system of first-order diff.equations.

If we let:

and write the above a system of ODEs:

The system will have the same number of equations as the order to the diff.equation.

We demonstrate by solving the "mathematical" pendulum (ie. ideal -> no losses -> perpetual motion!)

For a pendulum of mass m, length L from pivot, at angle ϕ, under influence of gravity g.

Starting angle

Starting angular velocity

Let

Finally, in vectorized form:

In the script, we use odeset() to give ode45() trigger events. In this way, we can find the length of a period, by terminating the solver once

(the angular velocity) falls to 0.

close all

f1 = @(u1, u2) u2; % The differential equations u' = f(t,u).

f2 = @(u1, u2, g, L) -g/L.*sin(u1);

f = @(t, u, g, L) [f1(u(1),u(2)); f2(u(1), u(2), g, L)]; % System of ode's

g = 9.81;

L = 0.1;

tspan = linspace(0, 2, 400);

phi0 = 30*pi/180;

u0 = [phi0; 0]; % Initial conditions.

%[t U] = ode45(@(t,u)f(t,u,g,L), tspan, u0);% solve ODE for u1 and u2. This gives us angle phi and angular velocity phi'.

opts = odeset('Events',@funevent_pendulum); % Event handler. Sets ode45() to react to certain events. See funevent() at bottom of script.

phi0_vec = (10:10:170).*pi./180; % A few initial values for starting angle phi0.

T_vec = zeros(size(phi0_vec));

i = 1;

for phi0 = phi0_vec

[t U] = ode45(@(t,u)f(t,u,g,L), tspan, [phi0;0], opts); % Starting conditions are angle and anglular velocity.

T_vec(i) = t(end);

i = i + 1;

subplot(121)

hold on

plot(t, U(:,1))

xlabel('t'); ylabel('φ(t)')

axis square

subplot(122)

plot(U(:,1), U(:,2))

hold on

xlabel('φ(t)'); ylabel('ω(t)')

axis square

end

figure

plot(phi0_vec, T_vec)

grid on

xlabel('φ_0');

ylabel('T [s]')

title('Pendulum period T vs starting angle φ_0')

Uppgift 6: Damped pendulum

The order of the diff.equation has not changed -> we still use a system of two equations.

We simply modify the expression for

close all

L = 0.1; m = 0.1; c = 0.2; g = 9.81;

f1 = @(u1, u2) u2; % The differential equations u' = f(t,u).

f2 = @(u1, u2, g, L, m) -g/L.*sin(u1) - c/m*u2;

f = @(t, u, g, L, m) [f1(u(1),u(2)); f2(u(1), u(2), g, L, m)]; % System of ode's

phi0 = deg2rad(50);

tspan = [0 1];

opts = odeset('Events',@funevent_damped_pendulum); % "Interrupt handler" for the ode solver. ie. quit when conditions met.

for phi0 = deg2rad(10:20:70)

[t, U] = ode45(@(t,u)f(t,u,g,L,m), tspan, [phi0;0]);

subplot(121)

plot(t, U(:,1))

hold on

subplot(122)

plot(U(:,1), U(:,2))

hold on

grid on

end

Randvärdesproblem/Boundary-value problems

Here we solve 1D finite-differences boundary-value problems of the type:

This is different to differential equations with only an initial value since we now have an end-bound. We solve this by forming a system of linear equations Au=b that solve u for each x.

First, we divide interval [a,b] uniformly into n+1 points (n intervals). ie. we discretize the problem space.

We approximate the 2nd derivative on the l.h.s. with a 2nd order central difference

. We approximate the 1st derivative on the r.h.s with an "integer spanning" 1st order central difference

. This results in:

As an example, we consider a stationary (time-invariant) heat conductor of length L along the x-axis. The temperature at each point depends on the endpoints (the boundary values).

We replace the derivatives with central-difference approximations and restrict ourselves to the "inner grid" which excludes the endpoints (where the central difference is undefined).

Clearly, we can now write an expression for each

on the inner grid.

This results in the linear system of equations:

Matrix A is sparse, in this case a tridiagonal matrix.

We can generate it using spdiags() (sparse diagonals). Think of it as a rectangular matrix that we then shear.

||| \\\

||| => \\\

||| \\\

close all

% Problem parameters

g0 = 40;

gL = 60;

L = 1;

N = 30;

k = 2;

f = @(x) 200*exp(-(x-L/2).^2);

% Specify problem space

h = L/(N+1);

i = 0:N+1;

x = i.*h;

x_inner = x(2:end-1);

% Generate tridiagonal A coefficient matrix.

ett = ones(N,1);

A = spdiags([-ett 2*ett -ett], [-1 0 1], N, N);

% ^ for row i, place the input pattern at col i-1, i, i+1

full(A); % Output full matrix (not in sparse form).

% Form b matrix, selecting only inner grid.

b = (h^2/k.*f(x_inner))';

% Insert boundary values. According to 2nd order central difference, the

% boundary values appear as elements u_0 and u_N+1 in the expressions for f_1 and f_N, respectively. We simply move them to the r.h.s. ie. matrix b.

b(1) = b(1) + g0;

b(end) = b(end) + gL;

% Solve for u on the inner grid.

u = A\b;

% Finally complete the solution with the boundary values

u = [g0; u; gL ];

plot(x, u)

title('-ku''''=200exp(-(x-L/2)), g_0=40, g_L=60')

xlabel('x')

ylabel('Temp(x)')

grid on

axis([0 1 40 65])

3.2 Inskjutningsmetoden/Shooting method (https://en.wikipedia.org/wiki/Shooting_method)

In the previous boundary-value problem, we rewrote a 2nd order differential equation in terms of Euler approximations and formed a large system of equations (one for each discrete point in the problem space [excluding end-points]).

With this method, we reduce the boundary-value problem to an initial-value problem. We vary the initial parameters ("shooting" along different trajectories) until we find a solution that satisfies the boundary conditions. Basically, it makes systemic the process of estimating initial values, solving and adjusting until a solution is found.

Consider the boundary value problem from before:

We can rewrite this as a system of 1st order equations:

Where g() is an arbitrary function that gives 0 when the boundary conditions are satisfied.

An equivalent initial-value problem would be:

If we write the solution as a function of x and the initial value

. That is, this is the function that will evolve as a result of the initial value.

The solution to the boundary-value problem is then:

In other words, for initial value s, the function will evolve in such a way that

is the ending boundary, thus solving the problem.

Numerically, we solve this by finding the root of q(s) eg. using fzero/newton's method.

Solve the previous heat transfer problem using the Shooting method

A note on functions, anonymous functions and function-handles

Anonymous function are limited to a single executable statement, beyond that we must use "real" functions.

If "real" functions call other user defined functions, anonymous functions won't be in scope (unless passed with a function-handle), and so they must also be "real" functions.

The '@' creates a function-handle. Anonymous functions are already defined as function-handles, but "real" functions are not, so they must be given one before they can be passed as an argument.

"Real" functions are special in that they can be called directly from anywhere without the need for a function handle.

close all

g0 = 40;

gL = 60;

L = 1;

k = 2;

s = linspace(-10,110); % Guesses for unknown u' that will result in correct upper boundary value u(b,s)=gL.

xspan = [0 1]; % Problem domain, NOT boundary conditions.

q = zeros(size(s));

for k = 1:length(s)

[x, u] = ode45(@ode_system, xspan, [g0; s(k)]);

q(k) = u(end,1) - gL; % Compare obtained boundary value with actual one.

end

plot(s,q)

hold on

title('Error in upper bound (q(s)) for various guesses for initial conditions (s) ')

xlabel('s')

ylabel('q(s)')

grid on

q_root = fzero(@(s)q_heat_transfer(s), 0) % Find initial value s that results in boundary satisfied.

q_root = 66.1281

plot(q_root, 0, 'ro')

figure

[x u] = ode45(@ode_system, xspan, [g0; q_root]);

plot(x, u(:,1))

title('Heat transfer')

xlabel('x'); ylabel('u(x)')

hold on

%plot(x, u(:,2))

grid on

axis([0 1 40 65])

function ret = f_source(x) % functions for uppgift 8.
    ret = 200*exp(-(x-1/2).^2);
end
function q = q_heat_transfer(s)
    [x u] = ode45(@ode_system, [0 1], [40; s]);
    q = u(end,1) - 60;
end
function uprim = ode_system(x,U)
    uprim = [U(2); f_source(x)./(-2)];
end

function q = findq(s) % functions uppgift 9
[t u] = ode45(@odesystem_uppg9, [0 1], [0; s]);
q = u(end,1) - 0;
end
function uprim = odesystem_uppg9(t,U) 
    uprim = [U(2); 1-U(1).^2];
end

function [value, isterminal, direction]=funevent_pendulum(t,U)
    value = U(2); % This is the variable to react to. Angular velocity. 
    isterminal = 1; % Quit ode45() once this event occurs.
    direction = -1; % +-1: ignore unless rising/falling towards zero. 0 to allow any. We want to angular velocity to rise towards zero (it swings back in positive direction to complete the period). 
end
function [value, isterminal, direction]=funevent_damped_pendulum(t,U)
i = 4
    value = U(2); % This is the variable to react to. Angular velocity. 
    isterminal = 1; % Quit ode45() once this event occurs.
    direction = -1; % +-1: ignore unless rising/falling towards zero. 0 to allow any. We want to angular velocity to rise towards zero (it swings back in positive direction to complete the period). 
end