How can I use Lorenz Attractor code?

140 次查看(过去 30 天)
Hi everyone!
i want to simulate Lorenz Attractor using the script I found in Matlab File Exchange by Moiseev Igor. But I do not know how to input my parametes here.
my parameters are
sigma= 10
beta=8/3;
rho=28
x=5
y=5
z=5
and i want to integrate it from t=0 to 20.
function [x,y,z] = lorenz(rho, sigma, beta, initV, T, eps)
% LORENZ Function generates the lorenz attractor of the prescribed values
% of parameters rho, sigma, beta
%
% [X,Y,Z] = LORENZ(RHO,SIGMA,BETA,INITV,T,EPS)
% X, Y, Z - output vectors of the strange attactor trajectories
% RHO - Rayleigh number
% SIGMA - Prandtl number
% BETA - parameter
% INITV - initial point
% T - time interval
% EPS - ode solver precision
%
% Example.
% [X Y Z] = lorenz(28, 10, 8/3);
% plot3(X,Y,Z);
if nargin<3
error('MATLAB:lorenz:NotEnoughInputs','Not enough input arguments.');
end
if nargin<4
eps = 0.000001;
T = [0 25];
initV = [0 1 1.05];
end
options = odeset('RelTol',eps,'AbsTol',[eps eps eps/10]);
[T,X] = ode45(@(T,X) F(T, X, sigma, rho, beta), T, initV, options);
plot3(X(:,1),X(:,2),X(:,3));
axis equal;
grid;
title('Lorenz attractor');
xlabel('X'); ylabel('Y'); zlabel('Z');
x = X(:,1);
y = X(:,2);
z = X(:,3);
return
end
function dx = F(T, X, sigma, rho, beta)
% Evaluates the right hand side of the Lorenz system
% x' = sigma*(y-x)
% y' = x*(rho - z) - y
% z' = x*y - beta*z
% typical values: rho = 28; sigma = 10; beta = 8/3;
dx = zeros(3,1);
dx(1) = sigma*(X(2) - X(1));
dx(2) = X(1)*(rho - X(3)) - X(2);
dx(3) = X(1)*X(2) - beta*X(3);
return
end
Thank you and have a nice day.
  2 个评论
KALYAN ACHARJYA
KALYAN ACHARJYA 2019-5-25
编辑:KALYAN ACHARJYA 2019-5-25
Is these x,y,z are same as function output [x,y,z]?
x=5
y=5
z=5
Function
function [x,y,z] = lorenz(rho, sigma, beta, initV, T, eps)
%.........^^^^
end
Darwin Tuazon
Darwin Tuazon 2019-5-25
I'm confused to function [x,y,z] = lorenz(rho, sigma, beta, initV, T, eps)
Parameters values are sigma= 10; beta=8/3; rho=28;
and my intial values is x=5; y=5; z=5

请先登录,再进行评论。

回答(3 个)

Sulaymon Eshkabilov
编辑:Sulaymon Eshkabilov 2019-5-26
Hi,
You were not executing the codes properly. Here is a single code that associates both scripts into one. Now it is much simpler.
sigma=10; beta=8/3; ro=28; % Your data
ICs=[5, 5, 5]; % Your data
t=[0, 20];
OPTs = odeset('reltol', 1e-6, 'abstol', 1e-8);
[time, fOUT]=ode45(@(t, x)([-sigma*x(1)+sigma*x(2); -x(2)-x(1).*x(3); -beta*x(3)+x(1).*x(2)-beta*ro]), t, ICs, OPTs);
close all
figure
plot3(fOUT(:,1), fOUT(:,2), fOUT(:,3)), grid
xlabel('x(t)'), ylabel('y(t)'), zlabel('z(t)')
title('LORENZ functions x(t) vs. y(t) vs. z(t)')
axis tight
figure
comet3(fOUT(:,1), fOUT(:,2), fOUT(:,3))
figure
subplot(311)
plot(time, fOUT(:,1), 'b','linewidth', 3), grid minor
title 'LORENZ functions x(t), y(t), z(t)', xlabel 'time', ylabel 'x(t)'
subplot(312)
plot( time', fOUT(:,2), 'r', 'linewidth', 2 ), grid minor
xlabel 'time', ylabel 'y(t)'
subplot(313)
plot(time, fOUT(:,3),'k', 'linewidth', 2), grid minor, xlabel 'time', ylabel 'z(t)'
figure
plot(fOUT(:,1), fOUT(:,2), 'b', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'y(t)'
axis square
figure
plot(fOUT(:,1), fOUT(:,3), 'k', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'z(t)'
axis square
figure
plot(fOUT(:,2), fOUT(:,3), 'm', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('y(t)'), ylabel 'z(t)'
axis square
Good luck.
  3 个评论
Sulaymon Eshkabilov
Sulaymon Eshkabilov 2021-11-14
Post your question here in: www.mathworks.com
Ghofran Khaled
Ghofran Khaled 2021-11-15
I posted it before and didn't get an answer
The question is in this link
https://ww2.mathworks.cn/matlabcentral/answers/1586514-how-can-i-create-improved-lorenz-mapping-code?s_tid=srchtitle_how%20can%20i%20create%20improved%20lorenz%20code_1

请先登录,再进行评论。


Sulaymon Eshkabilov
Hi Darwin,
Here is my version of the Lorenz Atractor simulation code:
function df = LORENZ_sys_1ODE(~, x)
% HELP: Lorenz Functions
% dx/dt=-sigma*x+sigma*y;
% dy/dt=- y-x*z;
% dz/dt=-beta*z+x*y-beta*ro;
sigma=10; beta=8/3; ro=28;
% ICs: x(0)=5; y(0)=5; z(0)=5; % Your ICs
df=[-sigma*x(1)+sigma*x(2); ...
-x(2)-x(1).*x(3);...
-beta*x(3)+x(1).*x(2)-beta*ro];
end
Run this part to simulate the whole system
ICs=[5, 5, 5]; % Your data
t=[0, 20]; % Your simulation space
OPTs = odeset('reltol', 1e-6, 'abstol', 1e-8); % Optional ODE options set up
[time, fOUT]=ode45(@LORENZ_sys_1ODE, t, ICs, OPTs);
close all
figure
plot3(fOUT(:,1), fOUT(:,2), fOUT(:,3)), grid
xlabel('x(t)'), ylabel('y(t)'), zlabel('z(t)')
title('LORENZ functions x(t) vs. y(t) vs. z(t)')
axis tight
figure
comet3(fOUT(:,1), fOUT(:,2), fOUT(:,3))
figure
subplot(311)
plot(time, fOUT(:,1), 'b','linewidth', 3), grid minor
title 'LORENZ functions x(t), y(t), z(t)', xlabel 'time', ylabel 'x(t)'
subplot(312)
plot( time', fOUT(:,2), 'r', 'linewidth', 2 ), grid minor
xlabel 'time', ylabel 'y(t)'
subplot(313)
plot(time, fOUT(:,3),'k', 'linewidth', 2), grid minor, xlabel 'time', ylabel 'z(t)'
figure
plot(fOUT(:,1), fOUT(:,2), 'b', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'y(t)'
axis square
figure
plot(fOUT(:,1), fOUT(:,3), 'k', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'z(t)'
axis square
figure
plot(fOUT(:,2), fOUT(:,3), 'm', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('y(t)'), ylabel 'z(t)'
axis square
Good luck.
  1 个评论
Darwin Tuazon
Darwin Tuazon 2019-5-26
Thank you for your response. I compile your code Sulaymon Eshkabilov but there is an error.
>> LORENZ_sys_1ODE
Error: File: LORENZ_sys_1ODE.m Line: 7 Column: 11
Incorrect use of '=' operator. To assign a value to a variable, use '='. To compare values for
equality, use '=='.
Is there anything I did wrong?
Please see your code below.
function df = LORENZ_sys_1ODE(~, x)
% % HELP: Lorenz Functions
% % dx/dt=-sigma*x+sigma*y;
% % dy/dt=- y-x*z;
% % dz/dt=-beta*z+x*y-beta*ro;
sigma=10; beta=8/3; ro=28;
ICs: x(0) =5; y(0) = 5; z(0) = 5; % Your ICs
df=[-sigma*x(1)+sigma*x(2); ...
-x(2)-x(1).*x(3);...
-beta*x(3)+x(1).*x(2)-beta*ro];
ICs=[5, 5, 5]; % Your data
t=[0, 20]; % Your simulation space
OPTs = odeset('reltol', 1e-6, 'abstol', 1e-8); % Optional ODE options set up
[time, fOUT]=ode45(@LORENZ_sys_1ODE, t, ICs, OPTs);
close all
figure
plot3(fOUT(:,1), fOUT(:,2), fOUT(:,3)), grid
xlabel('x(t)'), ylabel('y(t)'), zlabel('z(t)')
title('LORENZ functions x(t) vs. y(t) vs. z(t)')
axis tight
figure
comet3(fOUT(:,1), fOUT(:,2), fOUT(:,3))
figure
subplot(311)
plot(time, fOUT(:,1), 'b','linewidth', 3), grid minor
title 'LORENZ functions x(t), y(t), z(t)', xlabel 'time', ylabel 'x(t)'
subplot(312)
plot( time', fOUT(:,2), 'r', 'linewidth', 2 ), grid minor
xlabel 'time', ylabel 'y(t)'
subplot(313)
plot(time, fOUT(:,3),'k', 'linewidth', 2), grid minor, xlabel 'time', ylabel 'z(t)'
figure
plot(fOUT(:,1), fOUT(:,2), 'b', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'y(t)'
axis square
figure
plot(fOUT(:,1), fOUT(:,3), 'k', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'z(t)'
axis square
figure
plot(fOUT(:,2), fOUT(:,3), 'm', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('y(t)'), ylabel 'z(t)'
axis square
end

请先登录,再进行评论。


Sulaymon Eshkabilov
Sulaymon Eshkabilov 2021-11-14
You simulated the code incorrectly. Here is how you should run the code in one m-file:
ICs=[5, 5, 5]; % Your data
t=[0, 20]; % Your simulation space
OPTs = odeset('reltol', 1e-6, 'abstol', 1e-8); % Optional ODE options set up
[time, fOUT]=ode45(@LORENZ_sys_1ODE, t, ICs, OPTs);
close all
figure
plot3(fOUT(:,1), fOUT(:,2), fOUT(:,3)), grid
xlabel('x(t)'), ylabel('y(t)'), zlabel('z(t)')
title('LORENZ functions x(t) vs. y(t) vs. z(t)')
axis tight
figure
comet3(fOUT(:,1), fOUT(:,2), fOUT(:,3))
figure
subplot(311)
plot(time, fOUT(:,1), 'b','linewidth', 3), grid minor
title 'LORENZ functions x(t), y(t), z(t)', xlabel 'time', ylabel 'x(t)'
subplot(312)
plot( time', fOUT(:,2), 'r', 'linewidth', 2 ), grid minor
xlabel 'time', ylabel 'y(t)'
subplot(313)
plot(time, fOUT(:,3),'k', 'linewidth', 2), grid minor, xlabel 'time', ylabel 'z(t)'
figure
plot(fOUT(:,1), fOUT(:,2), 'b', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'y(t)'
axis square
figure
plot(fOUT(:,1), fOUT(:,3), 'k', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'z(t)'
axis square
figure
plot(fOUT(:,2), fOUT(:,3), 'm', 'linewidth', 1.5)
grid minor, title('LORENZ functions'), xlabel('y(t)'), ylabel 'z(t)'
axis square
function df = LORENZ_sys_1ODE(~, x)
% HELP: Lorenz Functions
% dx/dt=-sigma*x+sigma*y;
% dy/dt=- y-x*z;
% dz/dt=-beta*z+x*y-beta*ro;
sigma=10; beta=8/3; ro=28;
% ICs: x(0)=5; y(0)=5; z(0)=5; % Your ICs
df=[-sigma*x(1)+sigma*x(2); ...
-x(2)-x(1).*x(3);...
-beta*x(3)+x(1).*x(2)-beta*ro];
end

类别

Help CenterFile Exchange 中查找有关 General Physics 的更多信息

产品

Community Treasure Hunt

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

Start Hunting!

Translated by