Description

Exercise 1: Initial Value Problems
Prof. Dr. Moritz Diehl, Florian Messerer, Andrea Zanelli, Dimitris Kouzoupis

In this exercise we will have a first look at initial value problems (IVPs). We will analyze some of the conditions necessary for the existence and uniqueness of solutions to IVPs and we will implement and compare explicit integration methods for the solution of ordinary differential equations.
1. Existence and uniqueness of solutions to IVPs: Lipschitz-continuity. In order to guarantee that a solution to an IVP exists and it is unique, we can rely on Theorem 1.2. Among the conditions that have to be satisfied for the results of the theorem to hold, the function f(x(t),t) describing the dynamics of the system, has to be Lipschitz continuous with respect to x.
(a) Indicate (without proof) which of the following functions are globally Lipschitz (GL), locally Lipschitz (LL) or not locally Lipschitz (NLL):
(GL) (LL) (NLL)
i. x2
ii. iii.
iv.
v. kxk2
vi.
vii. kAxk2
viii. sinxTAx
ix. sinkxk2
2. Explicit integrators. In Section 1.2 we will introduce integration schemes for ordinary differential equations (ODEs). These schemes allow one to compute a numerical approximation of the true solution to a initial value problem (IVP) of the form
x˙(t) = f(x(t)), t ∈ [0,T]
(1) x(0)= x0.
In this exercise we will compare three different explicit integration schemes.
(a) The MATLAB function ode45 provides a user-friendly interface to an implementation of the Runge-Kutta integrator of type Dormand-Prince. It is suitable for the integration of non-stiff ODEs of the form ˙x = f(t,x) and it can provide reasonable accuracy in reasonable computation times in many practical cases. Write a script that calls ode45 in order to obtain a numerical approximation of the solution to the IVP for a damped oscillator
x˙1(t) = x2(t)
x˙2(t) x(0) = −0.2×2 − x1(t)
= x0 (2)
1
for x0 = [1, 0]T and provides output on the grid defined by t = kTs, where k = 0,1,…,19 and Ts = 0.5. Plot the obtained trajectories.
(2 points)
(b) We will now implement two other explicit integration schemes. Firstly, implement the explicit Euler scheme
xn+1 = xn + hf(xn). (3)
Secondly, implement the explicit Runge-Kutta scheme of order four (RK4):
k1 = f(xn)
(4)
Using the two schemes (h = 0.5 for the RK4 scheme and h = 0.125 for explicit Euler), compute numerical approximations of the solution to the IVP in (2) and plot the obtained trajectories in the same figure as (a).
(c) Implement an RK4 integrator for our IVP as a CasADi function and use it to perform the numerical simulation.
Hint: you can build an expression from CasADi variables and then use it to declare a function, as in the following example. Use full() to convert a result to a normal Matlab variable.
1 import casadi.*
2 x = MX.sym(‘x’,2,1);
3 expr = sin(x(1))*x(2);
4 f = Function(‘f’, {x}, {expr});
5 y = f([pi/2;2]);
6 y2 = full(f([pi/2;2]));
3. The Lorenz attractor. We will now use the integrators implemented to simulate some interesting systems.
(a) In 1963, Edward Lorenz developed a simplified mathematical model for atmospheric convection. The model is a system of three ordinary differential equations now known as the Lorenz equations:
x˙1 = σ(x2 − x1)
x˙2 x˙3 = x1(ρ − x3) − x2
= x1x2 − βx3 (5)
with ρ = 28, σ = 10 and .
i. Does the Lorenz attractor satisfy the conditions necessary for the Picard-Lindeloef theorem to be applicable?
ii. Using the RK4 integrator implemented before, simulate the system for x0 = [1, 0, 0]T and t ∈ [0,100] with step-size h = 0.01. Plot the resulting trajectory in 3D (plot3()).
2