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

27 vues (au cours des 30 derniers jours)
Ughur
Ughur le 24 Fév 2025
Commenté : Voss le 24 Fév 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 commentaires
Walter Roberson
Walter Roberson le 24 Fév 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 le 24 Fév 2025
Thank you. I uploaded the code in pictures.

Connectez-vous pour commenter.

Réponse acceptée

Voss
Voss le 24 Fév 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;
  2 commentaires
Ughur
Ughur le 24 Fév 2025
Thank you so much !!
Voss
Voss le 24 Fév 2025
You're welcome!

Connectez-vous pour commenter.

Plus de réponses (0)

Catégories

En savoir plus sur Programming dans Help Center et File Exchange

Produits


Version

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by