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

16 visualizaciones (últimos 30 días)
Ughur
Ughur el 24 de Feb. de 2025
Comentada: Voss el 24 de Feb. de 2025
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
Walter Roberson
Walter Roberson el 24 de Feb. de 2025
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.
Ughur
Ughur el 24 de Feb. de 2025
Thank you. I uploaded the code in pictures.

Iniciar sesión para comentar.

Respuesta aceptada

Voss
Voss el 24 de Feb. de 2025
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 Help Center y File Exchange.

Productos


Versión

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by