Question

Implement the Lorenz-63 model in MATLAB, and solve numerically (using MATLAB’s ode45 or other built-in solver) for any random, non-zero initial conditions to reproduce, qualitatively, Figures 3 and 4. Note that, as usual, you should label your axes and make the plots as “pretty” as possible.

The model equations are

dx/ dt = σ(y − x),

dy /dt= x(ρ−z)−y,

dz/dt = xy − βz.
Use parameter values σ = 10, ρ = 28, and β = 8/3.

Answer 1

Define function in Matlab as below:

function lorenz

%clear all the pre defined variables and parameters

clear all

%clear screen

clc

%define the time scale

[0 1000];

%define the initial conditions for the lorenz system

x0 = [0.1, 0.1, 0.1];

%call the function as using ode45 method

[T,X] = ode45(@system,[0 100],x0);

%Plot the variables with respect to time

plot(X(:,1),'r','LineWidth',1.5)

hold on

plot(X(:,2),'g','LineWidth',1.5)

hold on

plot(X(:,3),'b','LineWidth',1.5)

%define title and labels of plot

xlabel('Time (t)','FontWeight','bold','fontsize',14) %x-axis label

ylabel('System','FontWeight','bold','fontsize',14) %y-axis label

title('Dynamic Behavior')%title

legend('x','y','z')

%Phase plot

figure

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

%define title and labels of plot

xlabel('x','FontWeight','bold','fontsize',14) %x-axis label

ylabel('y','FontWeight','bold','fontsize',14) %y-axis label

title('x-y Phase Plane')%title

figure

plot(X(:,1),X(:,3),'g','LineWidth',1.5)

%define title and labels of plot

xlabel('x','FontWeight','bold','fontsize',14) %x-axis label

ylabel('z','FontWeight','bold','fontsize',14) %y-axis label

title('x-z Phase Plane')%title

figure

plot(X(:,2),X(:,3),'b','LineWidth',1.5)

%define title and labels of plot

xlabel('y','FontWeight','bold','fontsize',14) %x-axis label

ylabel('z','FontWeight','bold','fontsize',14) %y-axis label

title('y-z Phase Plane')%title

%Plot 3D

figure

plot3(X(:,1),X(:,2),X(:,3),'LineWidth',1.5)

xlabel('x','FontWeight','bold','fontsize',14) %x-axis label

ylabel('y','FontWeight','bold','fontsize',14) %y-axis label

zlabel('z','FontWeight','bold','fontsize',14) %z-axis label

title('3D Lorenz Attractor')%title

%define the model function

%assume x = x1, y = x2 and z = x3

function dx = system(t,x)

dx = zeros(3,1);

%define parametes values for the system

sigma = 10; rho = 28; beta = 8/3;

%ODE system

dx(1) = sigma * x(2) - sigma * x(1);

dx(2) = rho * x(1) - x(1) * x(3) -x(2);

dx(3) = x(1) * x(2) - beta * x(3);

return