RK4 method to solve Lorenz attractor with error calculations O(k^4)

Question

Using RK4 solve the lorenz attractor system of ODEs.  Plot the error at t=2 as a function of timestep in loglog.
I believe I have the code correct for the RK4 to solve the ODE.  I get a solution that matches the expected value.  when I check my error it looks like the system is an O(k^1) forward euler type solution.  Clearly something is wrong but I cant find the mistake. 
clc;
clear all;

% initialize variables
% run time 0 to 100 divided into .01 (h) increments 
t(1) = 0;  
x(1) = 1;
y(1) = 1;
z(1) = 1;
c_e(1)= 1; % cauchy error
p(1) =1;
m=1;


% the ODE linked system to solve
% variable = @(linked variables)
%ODE in 3d is a function of x, y, z and t
a=@(t, x, y, z) 9*(y - x);  
b=@(t, x, y, z) x*(26-z) - y;
c=@(t, x, y, z) x*y -z;    %since beta is 1 I could simplify

%first loop changes the time step
for m = 1:5; 
%set step size
h= .01*(1/(2^m));
h2=h/2;
t= 0:h:100; 
%make the loop
for i = 1:(length(t)-1)    %run the loop from 1 to 99, using t so you dont have to dig down to change the length of time
    
    %solve for each of the coefficents
    % the first is straight forward
    %where is the point right now
    a1x = a(t(i), x(i), y(i), z(i));
    a1y = b(t(i), x(i), y(i), z(i));
    a1z = c(t(i), x(i), y(i), z(i));
    %solve for the second coefficent
    %move it forward a half step and average with the last position
    a2x = a(t(i), x(i)+(h2)*a1x, y(i)+(h2)*a1y, z(i)+(h2)*a1z); 
    a2y = b(t(i), x(i)+(h2)*a1x, y(i)+(h2)*a1y, z(i)+(h2)*a1z);
    a2z = c(t(i), x(i)+(h2)*a1x, y(i)+(h2)*a1y, z(i)+(h2)*a1z);
    %solve for the third coefficent
    %move it forward a half step and average with the last position
    a3x = a(t(i), x(i)+(h2)*a2x, y(i)+(h2)*a2y, z(i)+(h2)*a2z);
    a3y = b(t(i), x(i)+(h2)*a2x, y(i)+(h2)*a2y, z(i)+(h2)*a2z);
    a3z = c(t(i), x(i)+(h2)*a2x, y(i)+(h2)*a2y, z(i)+(h2)*a2z);
    %solve for the end point/4th coefficent
    %move it to end point of this step
    a4x = a(t(i), x(i)+(h*a3x),y(i)+(h*a3y), z(i)+(h*a3z));
    a4y = b(t(i), x(i)+(h*a3x),y(i)+(h*a3y), z(i)+(h*a3z));
    a4z = c(t(i), x(i)+(h*a3x),y(i)+(h*a3y), z(i)+(h*a3z));
    
    % put it all together with the weights for each of the coefficents
    x(i+1) = x(i)+(h/6)*(a1x+2*a2x+2*a3x+a4x);
    y(i+1) = y(i)+(h/6)*(a1y+2*a2y+2*a3y+a4y); 
    z(1+i) = z(i)+(h/6)*(a1z+2*a2z+2*a3z+a4z);
    
    
end

% f1 = figure ('name', 'Lorenz ODE');
%        plot3(x,y,z);
    q(:,m)=[abs(x(i)-x(i-1)) abs(y(i)-y(i-1)) abs(z(i)-z(i-1))];  
    %save error for a pair of adjacent points for later
    steps(:,m)=[h/(.01)];                                             
    %save the timestep value
   
end

for ii=1:m;
    c_e(:,ii)=[mean(q(:,ii))];                  
    %what is the average of the x,y,z error for each time step run
    %taken at the last time step
   
end
    
    
f2 = figure ('name', 'cauchy error');
   loglog(steps, c_e); grid on
 
f3 = figure ('name', 'cauchy error in x y z');
   loglog(steps, q); grid on

RK4 method to solve Lorenz attractor with error calculations O(k^4)

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

RK4 method to solve Lorenz attractor with error calculations O(k^4)

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论