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.

LorenzRK4

testnew51

%   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: 
% http://www.mathworks.com/matlabcentral/answers/94495-how-can-i-create-animated-gif-images-in-matlab
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)
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s