Something is wrong with the code. It must display a chaotic graph.
15 views (last 30 days)
Show older comments
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;
% 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 Comments
Walter Roberson
on 24 Feb 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.
Accepted Answer
Voss
on 24 Feb 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;
% Plotting
figure(2);
plot3(y,x,z)
xlabel('x'),ylabel('y'),zlabel('z'),legend(' ζ(s) = 0.97 + 0.03 cos(s/10)')
grid on;
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
