Exercise 6: Non-linear problems

Illustrative example

Consider the differential equation describing the motion of frictionless pendulum

where θ is the angle, g is acceleration of gravity and l is the length. Since θ depends only on t, (1) is an initial value problem. We can convert it into a system of first order ODEs and solve it numerically with the methods we have seen during the course. Introducing the angular velocity

the first order system looks like

Defining

we can discretize the problem with implicit Euler scheme as

and solve it with either fixed-point algorithm or Newton method.

Fixed-point algorithm

For better readability we define

and rearrange (3) to the fixed-point form as

The MatLab code below solves (3) with the fixed-point algorithm.

% Parameters

l=1;

g=9.81;

% Initial condition

U = [ pi/4 ; 0 ];

% Function f

f = @(U) [ U(2) ; -g/l*sin(U(1)) ];

% Options for time-stepping

dt = 0.01;

tmax = 4;

% Options for non-linear solver

itmax = 1000;

tol = 1e-5;

relaxation = 0.1;

% Initialization

t = 0;

V = U;

K = round((tmax-t)/dt);

theta = NaN(K,1);

omega = theta;

tic;

% Time-stepping loop

for k=1:K

% Fixed-point algorithm

error = inf;

it = 0;

while error>tol && itmax>it

Vold = V;

Vnew = U + dt*f(V);

V = relaxation*Vnew + (1-relaxation)*V;

it = it+1;

error = norm(V-Vold);

end

if it==itmax, disp('Not converged'); end

% Proceed to next time-step

U = V;

t = t + dt;

% Save trajectory

theta(k) = U(1);

omega(k) = U(2);

end

ComputationTime = toc

ComputationTime = 9.9517

% Plot phase diagram

plot(theta,omega);

Note how the implicit Euler scheme dissipates energy, although the original equation models a frictionless pendulum.

Newton-Raphson method

We can also re-arrange (3) into the form

and approximate the roots of the non-linear function

at each time step by Newton-Raphson method. See the MatLab script below:

% Re-initialization