Something is wrong with the code. It must display a chaotic graph.

clc; clear; close all;
% Parameters
h = 0.005; t(1)=0; tfinal = 10;
t = t(1):h:tfinal;
N = ceil((tfinal - t(1)) / h);
x(1) = 32;
y(1) = 32;
z(1) = 32;
alpha =0.995 ;
a = 10; m = -810;
k = 30; j = 8/3;
c = 28; l = 35.5;
r = 900;
% Functions
g1 = @(t, x, y, z) a.*(y - x);
g2 = @(t, x, y, z) c.*x - y + k.*z - x.*z - m;
g3 = @(t, x, y, z) -k.*x - l.*y - j.*z + x.*y + r;
% Main Loop
x(2) = x(1) + h * g1(t(1), x(1), y(1), z(1));
y(2) = y(1) + h * g2(t(1), x(1), y(1), z(1));
z(2) = z(1) + h * g3(t(1), x(1), y(1), z(1));
x(3) = x(2) + h * g1(t(2), x(2), y(2), z(2));
y(3) = y(2) + h * g2(t(2), x(2), y(2), z(2));
z(3) = z(2) + h * g3(t(2), x(2), y(2), z(2));
tic;
for p = 3:N
x(p+1) = x(p) + ((2 - alpha)/2).*((1 - alpha).*(g1(t(p), x(p), y(p), z(p)) - g1(t(p-1), x(p-1), y(p-1), z(p-1)))+alpha.*( ...
(23/12)*g1(t(p), x(p), y(p), z(p)).*h - (4/3)*g1(t(p-1), x(p-1), y(p-1), z(p-1)).*h ...
+ (5/12)*g1(t(p-2), x(p-2), y(p-2), z(p-2))));
y(p+1) = y(p) + ((2 - alpha)/2).*((1 - alpha).*(g2(t(p), x(p), y(p), z(p)) - g2(t(p-1), x(p-1), y(p-1), z(p-1)))+alpha.*( ...
(23/12)*g2(t(p), x(p), y(p), z(p)).*h ...
- (4/3)*g2(t(p-1), x(p-1), y(p-1), z(p-1)).*h ...
+ (5/12)*g2(t(p-2), x(p-2), y(p-2), z(p-2))));
z(p+1) = z(p) + ((2 - alpha)/2).*((1 - alpha).*(g3(t(p), x(p), y(p), z(p)) - g3(t(p-1), x(p-1), y(p-1), z(p-1)))+alpha.*( ...
(23/12)*g3(t(p), x(p), y(p), z(p)).*h ...
- (4/3)*g3(t(p-1), x(p-1), y(p-1), z(p-1)).*h ...
+ (5/12)*g3(t(p-2), x(p-2), y(p-2), z(p-2))));
t(p+1)=t(p)+h;
end
toc;
Elapsed time is 0.016362 seconds.
% Plotting
figure(2);
plot3(y,x,z)
xlabel('x'),ylabel('y'),zlabel('z'),legend(' ζ(s) = 0.97 + 0.03 cos(s/10)')
grid on;
In this code, I tried to graph a chaotic system using Newton interpolation. I think something is wrong with the 2nd and 3rd initial values. I tried to fix it using the Euler method, but it always produces the same linear graph.

2 comentarios

Please post the code itself instead of an image of the code. There are no released versions of MATLAB that are able to execute images of code.
Thank you. I uploaded the code in pictures.

Iniciar sesión para comentar.

 Respuesta aceptada

Here's a guess. The only change is putting ".*h" in a few places where it seemed likely to be missing.
clc; clear; close all;
% Parameters
h = 0.005; t(1)=0; tfinal = 10;
t = t(1):h:tfinal;
N = ceil((tfinal - t(1)) / h);
x(1) = 32;
y(1) = 32;
z(1) = 32;
alpha =0.995 ;
a = 10; m = -810;
k = 30; j = 8/3;
c = 28; l = 35.5;
r = 900;
% Functions
g1 = @(t, x, y, z) a.*(y - x);
g2 = @(t, x, y, z) c.*x - y + k.*z - x.*z - m;
g3 = @(t, x, y, z) -k.*x - l.*y - j.*z + x.*y + r;
% Main Loop
x(2) = x(1) + h * g1(t(1), x(1), y(1), z(1));
y(2) = y(1) + h * g2(t(1), x(1), y(1), z(1));
z(2) = z(1) + h * g3(t(1), x(1), y(1), z(1));
x(3) = x(2) + h * g1(t(2), x(2), y(2), z(2));
y(3) = y(2) + h * g2(t(2), x(2), y(2), z(2));
z(3) = z(2) + h * g3(t(2), x(2), y(2), z(2));
tic;
for p = 3:N
x(p+1) = x(p) + ((2 - alpha)/2).*((1 - alpha).*(g1(t(p), x(p), y(p), z(p)) - g1(t(p-1), x(p-1), y(p-1), z(p-1)))+alpha.*( ...
(23/12)*g1(t(p), x(p), y(p), z(p)).*h - (4/3)*g1(t(p-1), x(p-1), y(p-1), z(p-1)).*h ...
+ (5/12)*g1(t(p-2), x(p-2), y(p-2), z(p-2)).*h));
% ^^^ I added this
y(p+1) = y(p) + ((2 - alpha)/2).*((1 - alpha).*(g2(t(p), x(p), y(p), z(p)) - g2(t(p-1), x(p-1), y(p-1), z(p-1)))+alpha.*( ...
(23/12)*g2(t(p), x(p), y(p), z(p)).*h ...
- (4/3)*g2(t(p-1), x(p-1), y(p-1), z(p-1)).*h ...
+ (5/12)*g2(t(p-2), x(p-2), y(p-2), z(p-2)).*h));
% ^^^ and this
z(p+1) = z(p) + ((2 - alpha)/2).*((1 - alpha).*(g3(t(p), x(p), y(p), z(p)) - g3(t(p-1), x(p-1), y(p-1), z(p-1)))+alpha.*( ...
(23/12)*g3(t(p), x(p), y(p), z(p)).*h ...
- (4/3)*g3(t(p-1), x(p-1), y(p-1), z(p-1)).*h ...
+ (5/12)*g3(t(p-2), x(p-2), y(p-2), z(p-2)).*h));
% ^^^ and this
t(p+1)=t(p)+h;
end
toc;
Elapsed time is 0.016403 seconds.
% Plotting
figure(2);
plot3(y,x,z)
xlabel('x'),ylabel('y'),zlabel('z'),legend(' ζ(s) = 0.97 + 0.03 cos(s/10)')
grid on;

Más respuestas (0)

Categorías

Más información sobre Programming en Centro de ayuda y File Exchange.

Productos

Versión

R2024a

Preguntada:

el 24 de Feb. de 2025

Comentada:

el 24 de Feb. de 2025

Community Treasure Hunt

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

Start Hunting!

Translated by