Discretize model dx dt = σx − σy, dy dt = xρ − xz − y, dz dt = xy − βz, using Runge-Kutta methods (of order 2 or order 4).

Using the code :

clear all

close all

clc

Beta = [10;28;8/3];

x0 = [10;10;10];% I.C

dt = 0.001; % time step size

tspan = 0:dt:20;

figure(1)

options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,3));

[t,x] = ode45(@(t,x)lorenz(t,x,Beta),tspan,x0,options);

plot3(x(:,1),x(:,2),x(:,3),'r','lineWidth',1.5);

set(gca,'color','w','xcolor','r','ycolor','r','zcolor','r')

set(gcf,'color','w')

xlabel('porportional to the rate of convention')

ylabel('horizontal temperature variation')

zlabel('vertical temperature variation')

hold on

options1 = odeset('RelTol',1e-6,'AbsTol',1e-6*ones(1,3));

[t1,x1] = ode45(@(t1,x1)lorenz(t1,x1,Beta),tspan,x0,options1);

plot3(x1(:,1),x1(:,2),x1(:,3),'b','lineWidth',1.5);

set(gca,'color','w','xcolor','b','ycolor','b','zcolor','b')

set(gcf,'color','w')

grid on

legend('1e-12','1e-6')

figure(2)

subplot(3,1,1)

plot(t,x(:,1),'r',t1,x1(:,1),'b')

title('x(t)') %porportional to the rate of convention in respect to time

legend('1e-12','1e-6')

xlabel('time')

ylabel('y(t)') %porportional to the rate of convention

subplot(3,1,2)

plot(t,x(:,2),'r',t1,x1(:,2),'b')

title('y(t)') %horizontal temperature variation in respect to time

ylabel('HA temp var')

subplot(3,1,3)

plot(t,x(:,3),'r',t1,x1(:,3),'b')

title('z(t)') %vertical temperature variation in respect to time

xlabel('time')

ylabel('VA temp var') %vertical temperature variation

find the Runge-Kutta methods (of order 2 or order 4) Use the same parameter values 10;28;8/3 and initial condition x(0) = 10, y(0) = 10, z(0) = 10, run the simulation for 0 ≤ t ≤ 20. Plot the approximated solution using Runge-Kutta methods and the solution using ode45 in the same figure.

[Walter Roberson]

What question are you asking of us? Are you getting an error?

Answer 1

[Karim]

After having a look at the wikipedia article (Lorenz system - Wikipedia), it appears that you have some mistakes in your equations that you showed in the title. The equations should be:

• dx/dt = σ *(y-x)
• dy/dt = x*(ρ−z) − y
• dz/dt = x*y − β*z

Based on your script Beta = [10; 28; 8/3], i assume that 'Beta' holds the values for sigma rho and beta? So

• σ = 10
• ρ = 28
• β = 8/3

the part missing from your code is the lorenz function itself. So we need to create a function. An example can be:

% create the anonymous Function
Beta = [10;28;8/3];
lorenz = @(t,x) [   Beta(1)*( x(2)-x(1) );          % this represents dx/dt
                    x(1)*( Beta(2)-x(3) ) - x(2);   % this represents dy/dt
                    x(1)*x(2) - Beta(3)*x(3)   ];   % this represents dz/dt

Now we can solve using the ode45 solver

% initial values
x0 =  [1;1;1];
% time span
dt = 1e-3;
tspan = 0:dt:20;

% set the options for the solver
options = odeset('RelTol',1e-12,'AbsTol',1e-12*ones(1,3));
% call the solver using the 'lorenz' function
[t,x] = ode45(lorenz,tspan,x0,options);

% create a figure
figure
plot3(x(:,1),x(:,2),x(:,3),'r','lineWidth',1.5);
xlabel('porportional to the rate of convention')
ylabel('horizontal temperature variation')
zlabel('vertical temperature variation')
grid on