Processing math: 100%
Francesco Tudisco

Associate Professor (Reader) in Machine Learning

School of Mathematics, The University of Edinburgh
The Maxwell Institute for Mathematical Sciences
Alan Turing Institute | Bayes Centre | Gen AI Lab School of Mathematics, Gran Sasso Science Institute email: f dot tudisco at ed.ac.uk

Lorenz attractor


Consider the following Lorenz system of ODEs y′(x)=[y′1y′2y3]=[−αy1+αy2,−y1y3+βy1−y2,y1y2−γy3.]=[f1(x,y1,y2,y3)f2(x,y1,y2,y3)f3(x,y1,y2,y3)]=f(x,y)

  1. Create a matlab function 

    function f = LORENZ(x,y,var)
    (on a separate file) that implements the ODE (1), with var containing the parameters (α,β,γ).

  2. Implement the implicit Adams method with 6 steps as a solver for the system (1). Use a predictor-corrector strategy, within a specific matlab function 

    function [x,y] = IAM6(FUNC,x0,xN,y0,N,var)
    where y0 is a vector that contains sufficiently many starting values. Use the following explicit Adams method for the predictor phase: yn+1=yn+h1440(4277fn−7923fn−1+9982fn−2−7298fn−3+2877fn−4−475fn−5) and the following implicit one for the corrector phase yn+1=yn+h1440(475fn+1+1427fn−798fn−1+482fn−2−173fn−3+27fn−4).

  3. Solve (1) for x∈[0,,30] with the parameters α=10, β=28, γ=8/3, N=10000 and y0=(1,1,1). Use the function RK6impl(FUNC,x0,xN,y0,N,var) provided below to compute the 6 starting values for the IAM6 method. Plot the resulting function on the three planes y1=0, y2=0 and y3=0, on three different figures, as shown below.

..

function [t, y] = RK6IMPL(FUNC, x0, xN, y0, N, var)

% RK6IMPL Runge-Kutta-Algorithm
%    Executes the implicit Runge-Kutta-Algorithm of Kuntzmann
%    and Butcher to solve a system of differential equations
%    of first order defined in by 'FUNC'.
%    FUNC	: Ordinary differential equation
%    x0 	: starting x
%    xN 	: end x
%    y0  	: initial value
%    N  	: number of x intervalls
%    var	: Function arguments passed to the ode file
%
% Butcher-Scheme of the implemented Runge-Kutta-Method:
%   1/2-\sqrt{15}/10  |      5/36          2/9-\sqrt{15}/15 5/36-\sqrt{15}/30
%        1/2          | 5/36+\sqrt{15}/24      2/9          5/36-\sqrt{15}/24
%   1/2+\sqrt{15}/10  | 5/36+\sqrt{15}/30  2/9+\sqrt{15}/15      5/36
%  ---------------------------------------------------------------------------
%                     |      5/18              4/9               5/18
%

tol = 10e-6;
MAXSTEP = 500;

% Computation of the coefficients
alpha =[ 0.5-sqrt(15)/10; 0.5; 0.5+sqrt(15)/10];
beta = [ 5/36, 2/9-sqrt(15)/15, 5/36-sqrt(15)/30;...
				 5/36+sqrt(15)/24, 2/9, 5/36-sqrt(15)/24;...
			   5/36+sqrt(15)/30, 2/9+sqrt(15)/15, 5/36];
gamma = [5/18, 4/9, 5/18];

h = (tN-t0)/N;
t(1) = t0;
y(:,1) = y0;

for i=1:N
	t(i+1) = t0 + i*h;
	% Initialization of k and l
	knew = zeros(length(y0),3);
	kold = ones(length(y0),3);
	l = 0;
	% Computation of k
	while (norm([knew(:,1);knew(:,2);knew(:,3)]...
              -[kold(:,1);kold(:,2);kold(:,3)]) > tol) & (l < MAXSTEP)
		kold = knew;
		knew(:,1) = feval(FUNC, t(i)+h*alpha(1),...
                      y(:,i)+h*(beta(1,1)*kold(:,1)+beta(1,2)*kold(:,2)...
                                +beta(1,3)*kold(:,3)), var);
		knew(:,2) = feval(FUNC, t(i)+h*alpha(2),...
                      y(:,i)+h*(beta(2,1)*kold(:,1)+beta(2,2)*kold(:,2)...
                                +beta(2,3)*kold(:,3)), var);
		knew(:,3) = feval(FUNC, t(i)+h*alpha(3),...
                      y(:,i)+h*(beta(3,1) * kold(:,1)+beta(3,2)*kold(:,2)...
                                +beta(3,3)* kold(:,3)), var);
		l = l+1;
	end
	% compute y_{i+1}
	y(:,i+1) = y(:,i) + h * knew * gamma;
end