This is the 4d Lorenz equation. how to solve it on Matlab?

dx/dt=y+y^2-a*y*z;

dy/dt=-z^2+b*y*z-u;

dz/dt=x*y;

du/dt=-c*x;

[John D'Errico]

Hint: Read the help for ODE45.

[Nix Jr]

Here has the function with code. But I didn't understand what's the problem. If you don't mind, would you please help me to fix the problem?

...............................................................function..........................................

function df = test(~, x)

% HELP: Lorenz Functions

%dx/dt=y+y^2-a*y*z;

%dy/dt=-z^2+b*y*z-u;

%dz/dt=x*y;

%du/dt=-c*x;

u=10; a=8/3; b=7; c=.5;

% ICs: x(0)=5; y(0)=5; z(0)=5; % Your ICs

df=[x(2)+([x(2)]^2)-a*x(2).*x(3); ...

-([x(3)]^2)+b*x(2).*x(3)-x(4);...

x(1).*x(2)...

-c*x(1)];

end

...................................................................................main code...........................

ICs=[5, 5, 5]; % Your data

t=[0, 20]; % Your simulation space

OPTs = odeset('reltol', 1e-6, 'abstol', 1e-8); % Optional ODE options set up

[time, fOUT]=ode45(@test, t, ICs, OPTs);

close all

figure

plot3(fOUT(:,1), fOUT(:,2), fOUT(:,3), grid

xlabel('x(t)'), ylabel('y(t)'), zlabel('z(t)')

title('LORENZ functions x(t) vs. y(t) vs. z(t)')

axis tight

figure plot3(fOUT(:,1), fOUT(:,2), fOUT(:,3),

figure

subplot(311)

plot(time, fOUT(:,1), 'b','linewidth', 3), grid minor

title 'LORENZ functions x(t), y(t), z(t)', xlabel 'time', ylabel 'x(t)'

subplot(312)

plot( time', fOUT(:,2), 'r', 'linewidth', 2 ), grid minor

xlabel 'time', ylabel 'y(t)'

subplot(313)

plot(time, fOUT(:,3),'k', 'linewidth', 2), grid minor, xlabel 'time', ylabel 'z(t)'

figure

plot(fOUT(:,1), fOUT(:,2), 'b', 'linewidth', 1.5)

grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'y(t)'

axis square

figure

plot(fOUT(:,1), fOUT(:,3), 'k', 'linewidth', 1.5)

grid minor, title('LORENZ functions'), xlabel('x(t)'), ylabel 'z(t)'

axis square

figure

plot(fOUT(:,2), fOUT(:,3), 'm', 'linewidth', 1.5)

grid minor, title('LORENZ functions'), xlabel('y(t)'), ylabel 'z(t)'

axis square

Answer 1

[Yongzhen Mi]

Hi, Nix Jr. I think the error is that you defined four variables in your differential equations, but you provided only three initial values to the ode45 function. Therefore, ICs should be something like this [5, 5, 5, 5]' (column vector better).

In addition, you can transfer parameters into the target function by using [time, fOUT]=ode45(@(x)test(~,x,a,b,c), t, ICs, OPTs);. In this way, the parameters a, b, and c can be claimed in the main function once, instead of being defined again and again when the target function is called.

[Nix Jr]

Thanks for you suggestion but this process doesn't work. or might be I didn't get your process.

[Yongzhen Mi]

Take this as an example:

a = 5;

b = 20;

c = 1;

d = 0.1;

k = 0.1;

e = 20.6;

h = 1;

F = @(t,Y) [a*(Y(2)-Y(1))-e*Y(4);Y(1)*Y(3)-h*Y(2);b-Y(1)*Y(2)-c*Y(3);k*Y(2)-d*Y(4)];

CI = [3.2 8.5 3.5 2.0]';

T = linspace(0,100,3000);

[t,Y]=ode45(F,T,CI);

plot3(Y(:,1),Y(:,2),Y(:,3))

I drew some pictures and the attractors could be observed. Please try again and hopefully this will work.