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

15 views (last 30 days)
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 Comments
Walter Roberson
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.

Sign in to comment.

Accepted Answer

Voss
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;
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;

More Answers (0)

Categories

Find more on Programming in Help Center and File Exchange

Products


Release

R2024a

Community Treasure Hunt

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

Start Hunting!