# Finding and Plotting Lorenz Solution using MATLAB

I use MATLAB to solve the following Lorenz initial value problem:
$\begin{cases} x'=-10(x+y) \\ y'=-x(z+28)-y \\ z'=xy-\frac{8}{3}z \\ x(0)=y(0)=z(0)=5 \end{cases}$

I wrote a function, LorenzRK4IVP(), that takes the system of three differential equations as input and solves the system using the Runge-Kutta method with step size $h=.01$. I plot the strange attractor as well as use MATLAB to produce a GIF of the solution.

%   Inputs:
%      f1,f2,f3 = y'(t,y) as a string
%            y0 = initial condition
%         inter = interval
%             h = step size
function [a] = LorenzRK4IVP(f1,f2,f3,inter,y0,h)
t(1)=inter(1);
y(1,:)=y0;
n=(inter(2)-inter(1))/h;

% evaluate derivative
function z=ydot(f1,f2,f3,t,y)
g1=inline(f1,'y1','y2','y3','t');
g2=inline(f2,'y1','y2','y3','t');
g3=inline(f3,'y1','y2','y3','t');
z(1)=g1(y(1),y(2),y(3),t(1));
z(2)=g2(y(1),y(2),y(3),t(1));
z(3)=g3(y(1),y(2),y(3),t(1));
end

% evaluate Runge-Kutta 4 formula
function y=trap(f1,f2,f3,t,y,h)

s1=ydot(f1,f2,f3,t,y);
s2=ydot(f1,f2,f3,t+(h/2),y+(h/2)*s1);
s3=ydot(f1,f2,f3,t+(h/2),y+(h/2)*s2);
s4=ydot(f1,f2,f3,t+(h),y+(h)*s3);

y=y+(h/6)*(s1+2*s2+2*s3+s4);
end

% loop
for i=1:n
t(i+1)=t(i)+h;
y(i+1,:)=trap(f1,f2,f3,t(i),y(i,:),h);
end
a=[t',y(:,1),y(:,2),y(:,3)];
end

% call function to solve Lorenz equations
L=LorenzRK4IVP('-10*y1+10*y2','-y1*y3+28*y1-y2','y1*y2-(8/3)*y3',[0,50],[5,5,5],.01)

% plot Lorenz solutions
plot3(L(:,2),L(:,3),L(:,4))

% create GIF source:
figure(1)
filename = 'Lorenz.gif';
for n = 1:(length(L)/10)
plot3(L(1:10*n,2),L(1:10*n,3),L(1:10*n,4))
drawnow
frame = getframe(1);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
if n == 1;
imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
else
imwrite(imind,cm,filename,'gif','WriteMode','append');
end
end

L=LorenzRK4IVP('-10*y1+10*y2','-y1*y3+28*y1-y2','y1*y2-(8/3)*y3',[0,50],[5,5,5],.01)


## 2 thoughts on “Finding and Plotting Lorenz Solution using MATLAB”

1. Royi Avital says:

Hi,

Any chance you make MCMC Blog posts with MATLAB code and not only R code?

1. AO says:

Nope sorry. I don’t use MATLAB too often. But the syntax is quite similar to R, so if you understand the R code/concept it shouldn’t be too bad to write a MATLAB translation.