CDS - Lecture 14 ~ Matlab Script 'sensitivity'

• Source: moodle

function sensitivity()
 
    close all
    
     Parameters
    b = 8/3;
    sigma = 10;
    r = 24.8;
 
     Find a solution
    tspan = 0:deltat:tFinal;
    
    % Find a solution using the first initial condition
    [t, u_1] = ode45(@f, tspan, [x0_1; y0_1; z0_1], [], sigma, r, b);
    x_1 = u_1(:,1);
    y_1 = u_1(:,2);
    z_1 = u_1(:,3);
    
    % Find a solution using the second initial condition
    [t, u_2] = ode45(@f, tspan, [x0_2; y0_2; z0_2], [], sigma, r, b);
    x_2 = u_2(:,1);
    y_2 = u_2(:,2);
    z_2 = u_2(:,3);    
    
     Plot the solution in the phase space 
    figure
    hold on
    plot3(x_1,y_1,z_1, 'b')
    plot3(x_2,y_2,z_2, 'r')
    title('Sensitivity to i.c. in the phase space')
    
    legend('Trajectory 1', 'Trajectory 2')
    xlabel('x')
    ylabel('y')
    zlabel('z')    
    % 
end
 
function du = f(t, u, sigma, r, b)
    x = u(1);
    y = u(2);
    z = u(3);
    
    dx = sigma*(y-x);
    dy = -x*z+r*x-y;
    dz = x*y-b*z;    
    
    du = [dx; dy; dz];
end