06_initial_value_ode

Problem 2

Part A

The analytical solution to y' = -y is y = e^-t

Part B

dt=[0.5 .1 .001];
figure();
hold on
     x=[0:dt(2):3]';
     y_n=zeros(size(x));
    y_n(1)=1;
for i=2:length(x)
    dy=prob2_ode(x,y_n(i-1));
    y_n(i)=y_n(i-1)+ dt(2)*dy;
end  
plot(x,exp(-x),x,y_n)
legend('Analytical','Euler','Location','Northeast')

Part C

function y=Prob2ptC(ode,dt,y0,tlimits)
maxit=100;
t=[tlimits(1):dt:tlimits(2)];
y=zeros(length(t),1);
y(1)=y0;
    for i=2:length(t)
        dy=ode(t(i-1),y(i-1));
        y(i)=y(i-1)+dt*dy;
        dyc=(dy+ode(t(i),y(i)))/2;
        y(i)=y(i-1)+dt*dyc;
        nits=1;

        while(1)
            nits=nits+1;
            yold=y(i);
            dyc=(dy+ode(t(i),y(i)))/2;
            y(i)=y(i-1)+dt*dyc;
            ea=abs(yold-y(i))/y(i)*100;
            if ea<.0001|nits>maxit;break;end;
        end;
    end
end

Problem 3

Part A

The analytical solution is y = cos(3t)

Part B

dt=[.1 .001];
figure();
hold on
x=[0:dt(2):3]';
y_n=zeros(length(x),2);
y_n(1,:)=[1,0];
for i=2:length(x)
        dy=prob3_ode(x,y_n(i-1,:));
        y_n(i,:)=y_n(i-1,:)+ dt(2)*dy;
end
plot(x,cos(3*x),x,y_n(:,1));

Part C

function y = prob3ptC(prob3_ode, dt, y0,tlimits)
    t = [tlimits(1):dt:tlimits(2)];
    y = zeros(length(t),length(y0));
    y(1,:) = y0;
    for i = 2:length(t)
        dy = prob3_ode(t(i-1),y(i-1,:));
        y(i,:) = y(i-1,:)+dt*dy;
        dyc = (dy + prob3_ode(t(i),y(i,:)))/2;
        y(i,:) = y(i-1,:)+dt*dyc;
        n = 1;
        while (1)
            n = n+1;
            yold = y(i,:);
            dyc = (dy+prob3_ode(t(i),y(i,:)))/2;
            y(i,:) = y(i-1,:)+dt*dyc;
            ea = abs(y(i,:)- yold)./y(i,:)*100;
            if max(ea) < 0.0001 || n>100; break, end
        end
    end
end

Problem 4

Part A

The analytical solution to v'=g-cd/m v*v is

Part B

cd = 0.25;
m = 60;
g = 9.81;
dt = 0.1;
t = [0:dt:12]';
x = zeros(length(t),2);
x(1,:) = [1000, 0];
for i = 2:length(t)
    dx = prob4_ode(t,x(i-1,:));
    x(i,:) = x(i-1,:)+dt*dx;
end

Part C

function x = prob4ptC(prob4_ode, dt, x0,tlimits)
    t = [tlimits(1):dt:tlimits(2)];
    x = zeros(length(t),length(x0));
    x(1,:) = x0;
    for i = 2:length(t)
        dx = prob4_ode(t(i-1),x(i-1,:));
        x(i,:) = x(i-1,:)+dt*dx;
        dxc = (dx + prob4_ode(t(i),x(i,:)))/2;
        x(i,:) = x(i-1,:)+dt*dxc;
        n = 1;
        while (1)
            n = n+1;
            yold = x(i,:);
            dxc = (dx+prob4_ode(t(i),x(i,:)))/2;
            x(i,:) = x(i-1,:)+dt*dxc;
            ea = abs(x(i,:)- yold)./x(i,:)*100;
            if max(ea) < 0.0001 || n>100; break, end
        end
    end
end

Problem 5

set(0, 'defaultAxesFontSize', 16)
set(0,'defaultTextFontSize',14)
set(0,'defaultLineLineWidth',1)
dt=0.1;
    x=[0:dt:20]';
    y_n=zeros(length(x),4);
    y_n(1,:)=[10 0 0 2];
    for i=2:length(x)
            dy=phugoid_ode(x(i-1),y_n(i-1,:));
            y_n(i,:)=y_n(i-1,:)+dt*dy;
    end

dt=0.01;
    x=[0:dt:20]';
    y_n2=zeros(length(x),4);
    y_n2(1,:)=[10 0 0 2];
    for i=2:length(x)
            dy=phugoid_ode(x(i-1),y_n2(i-1,:));
            y_n2(i,:)=y_n2(i-1,:)+dt*dy;
    end
    [t23,y23]=ode23(@(t,y) phugoid_ode(t,y),[0 20],[10 0 0 2]);
    d = figure(4);
    plot(y_n(:,3),y_n(:,4),y_n2(:,3),y_n2(:,4),y23(:,3),y23(:,4))
    legend('Euler .1','Euler .01','ode23','Location','Northeast')
    xlabel('Distance')
    ylabel('Height')
    ```