Hi everyone,
I started a thread before called "Lorenz Equations MATLAB", and this is based on the question, but im having a different problem.
So from before,
I have the following system of differential equations
I need to calculate the trajectory of in MATLAB using Heun's method (which I think is also called the trapezoidal method) with step size h=0.001, and initial conditions for t < 20.
Here is my matlab code for carrying this out (contains 3 seperate mfiles)
mfile 1
----------------------------------------------------------------------
function z = ydot(t,y)
% Lorenz Equations
% z: derivative of state y
s=10; r=28; b=8/3;
z=zeros(3,1);
z(1) = s*(y(2) - y(1));
z(2) = y(1)*(r-y(3))-y(2);
z(3) = y(1)*y(2) - b*y(3);
----------------------------------------------------------------------
mfile 2
----------------------------------------------------------------------
function y = trapezoid(y,h)
z1 = ydot(0,y);
g=y+h*z1';
z2=ydot(0,g);
y=y+ 0.5*h*(z1+z2)';
---------------------------------------------------------------------
mfile 3
---------------------------------------------------------------------
function [t,y] = lorentz(int,ic,h)
%
% this function solves and plots the Lorenz equation
% t : time step vector
% y : state variable (coordinates)
%
% int: integration interval
% ic: initial state coordinates
% h: integration step size
t = 0:h:20; % set up integration time steps
n = numel(t);
y = zeros(n,3); % preallocate space
y(1, : )=ic;
t(1) = int(1);
for k=2:n
y(k, : ) = trapezoid(y(k-1, : ),h);
end
figure(1); %bring figure 1 to the front
subplot(1,2,1)
plot3(y(:,1),y(:,2 ),y(:,3));
xlabel('x');
zlabel('z');
ylabel('y');
subplot(1,2,2)
plot(y(:,1),y(:,2)),grid on;
xlabel('y1');ylabel('y2');
title('Lorenz Attractor');
-----------------------------------------------------------------------
Now to run this, in the command window type
lorentz([0 20], [5 5 5], 0.001)
What I have to do now is calculate the difference between the two trajectories and .
I thought I could do this by typing in the command window
lorentz([0 20], [5 5 5], 0.001) - lorentz([0 20], [(5 + 0.0001) 5 5], 0.001)
but this returns a value of zero.
Any ideas on how to do this? Thanks in advance.