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).

7 次查看(过去 30 天)
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.

回答(1 个)

Karim
Karim 2023-3-15
编辑:Karim 2023-3-15
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

类别

Help CenterFile Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by