06_initial_value_ode

Homework 6

Problem 2

%This function which is utilized for Problem 2 of the 6th homework dies
%calculate the solution to an ode using the Euler and Heun methods
%Analytic solution (Part a)
ya=@(t) -exp(-t);

%Euler Method (Part b)
y0=1;
dt=0.01;
t=0:dt:3;
y_euler=zeros(1,length(t));
y_euler(1)=y0;
for j=2:length(t)
    y_euler(j)=y_euler(j-1)+dt*(-y_euler(j-1));
end
plot(t,y_euler)
hold on

%Heun method (Part c) 
y_heun=zeros(size(t));
dy=zeros(size(t));
dyc=zeros(size(t));
y_heun(1)=1;
n=0;
error=1;
for i=1:length(t)-1
    dy(i)=ya(t(i));
    dyc(i)=ya(t(i));
    y_heun(i+1)=y_heun(i)+dyc(i)*(t(i+1)-t(i));
    while error< 0.001 || n>100
        n=n+1;
        y_prev=y_heun(i+1);
        dyc(i)=(dy(i)+ya(t(i+1),y_heun(i+1)))/2;
        y_heun(i+1)=y_heun(i)+dyc(i)*(t(i+1)-t(i));
        error=abs(y_heun(i+1)-y_prev)/y_heun(i+1)*100;
    end
end

plot(t,y_heun)

Problem 3

%This function which is utilized for Problem 3 of the 6th homework dies
%calculate the solution to an ode using the Euler and Heun methods
%Analytic solution (Part a)
ya=@(t) cos(3*t);

%Euler Method (Part b)
y0=1;
dy0=0;
dt=0.01;
t=0:dt:3;
y_euler=zeros(1,length(t));
dy_euler=zeros(1,length(t));
y_euler(1)=y0;
dy_euler(1)=dy0;
for i = 1:length(t)-1
    dy_euler(i+1)=dy_euler(i)-9*y_euler(i)*dt;
    y_euler(i+1)=y_euler(i)+dt*dy_euler(i);
end
plot(t,y_euler)
hold on

%Heun method (Part c)
y_heun=zeros(size(t));
dy=zeros(size(t));
dyc=zeros(size(t));
y_heun(1)=1;
n=0;
error=1;

for i=2:length(t)-1
    dy(i)=dy(i-1)-9*y_heun(i-1)*dt;
    y_heun(i)=y_heun(i-1)+dt*dy(i-1);
    dyc(i)=(dy(i-1)+dy(i))/2;
    y_heun(i)=y_heun(i-1)+dt*dy(i);
    y_heun(i+1)=y_heun(i)+dyc(i)*(t(i+1)-t(i));
    while (1)
        n=n+1;
        y_prev=y_heun(i+1);
        dyc(i)=(dy(i)+ya(t(i+1)))/2;
        y_heun(i+1)=y_heun(i)+dyc(i)*(t(i+1)-t(i));
        error=abs(y_heun(i+1)-y_prev)/y_heun(i+1)*100;
        if error< 0.0001 || n>100 
            break
        end
    end
end
plot(t,y_heun)

Problem 4

%This function which is utilized for Problem 4 of the 6th homework dies
%calculate the solution to an ode using the Euler and Heun methods
%Analytic solution (Part a)
x0=1000;
g=9.81;
cd=0.25;
m=60;
dt=0.01;
t=0:dt:12;
ya=@(t) x0-(m/cd)*log(cosh(-sqrt((cd/m)*g)*t));

%Euler method (Part b)
x_euler=zeros(1,length(t));
x0=1000;
x_euler(1)=x0;
dx_euler=zeros(1,length(t));
v0=0;
dx_euler(1)=v0;

for i=2:length(t)
    x_euler(i)=x_euler(i-1)+dt*dx_euler(i-1);
    dx_euler(i)=dx_euler(i-1)-dt*(g-(cd/m)*(dx_euler(i-1))^2);
end
plot(t,x_euler)
hold on

%Heun method (Part c)
x_heun=zeros(size(t));
dx=zeros(size(t));
dxc=zeros(size(t));
dx(1)=v0;
dxc(1)=v0;
x_heun(1)=x0;
n=0;
error=1;

for i=2:length(t)-1
    dx(i)=dx(i-1)-dt*(g-(cd/m)*(dx(i-1))^2);
    x_heun(i)=x_heun(i-1)+dt*dx(i-1);
    dxc(i)=(dx(i-1)+dx(i))/2;
    x_heun(i)=x_heun(i-1)+dt*dx(i);
    x_heun(i+1)=x_heun(i)+dxc(i)*(t(i+1)-t(i));
    while (1)
        n=n+1;
        y_prev=x_heun(i+1);
        dxc(i)=(dx(i)+ya(t(i+1)))/2;
        x_heun(i+1)=x_heun(i)+dxc(i)*(t(i+1)-t(i));
        error=abs(x_heun(i+1)-y_prev)/x_heun(i+1)*100;
        if error< 0.0001 || n>100 
            break
        end
    end
end
plot(t,x_heun)

Problem 5

%This is Homework Six Problem 5 which calculates the vekiocity, angle, x
%and y position of a paper airplane when given inputs about its weight and
%drag coefficient. These properties are calculat3ed using the Euler method
%and the buiilt in ODE23 function.
%Part a

DL=1/5.2;
vt=5.5;
dt=0.1;
t=0:dt:10;
g=9.81;
y=zeros(4,length(t));
y(1,1)=10;
y(2,1)=3;
y(3,1)=0;
y(4,1)=2;
for i=2:length(t)
    y(1,i)=y(1,i-1)+dt*(-g*sin(y(2,i-1))-DL*(g/((vt)^2))*(y(1,i-1))^2);
    y(2,i)=y(2,i-1)+dt*(-(g/y(1,i-1))*cos(y(2,i-1))+(g/((vt)^2))*y(1,i-1));
    y(3,i)=y(3,i-1)+dt*y(1,i-1)*cos(y(2,i-1));
    y(4,i)=y(4,i-1)+dt*y(1,i-1)*sin(y(2,i-1));
end
%call function
dy=phugoid_ode(t,y)
%Part b

DL=1/5.2;
vt=5.5;
dt=0.01;
t=0:dt:20
g=9.81;
v=zeros(1,length(t));
v(1)=10;
theta=zeros(1,length(t));
theta(1)=0;
x=zeros(1,length(t));
x(1)=0;
y=zeros(1,length(t));
y(1)=2;
for i=2:length(t)
    v(i)=v(i-1)+dt*(-g*sin(theta(i-1))-DL*(g/((vt)^2))*(v(i-1))^2);
    theta(i)=theta(i-1)+dt*(-(g/v(i-1))*cos(theta(i-1))+(g/((vt)^2))*v(i-1));
    x(i)=x(i-1)+dt*v(i-1)*cos(theta(i-1));
    y(i)=y(i-1)+dt*v(i-1)*sin(theta(i-1));
end

plot(x,y)


dv=@(v,theta) -g*sin(theta)-DL*(g/vt^2)*v^2;
dtheta=@(theta,v) -(g/v)*cos(theta)+(g/vt^2)*v;
dX=@(v,theta) v*cos(theta);
dY=@(v,theta) v*sin(theta);
[t,X]=ode23(dX,[0,12],1000);




%Extra Credit


[v0,theta0]=meshgrid(linspace(0,30,20),linspace(0,2*pi,20));

for i=length(v0)
    for j=1:length(theta0)
        v=zeros(length(v0),length(theta0),length(t));
        v(i,j,1)=v0(i,1);
        theta=zeros(length(v0),length(theta0),length(t));
        theta(i,j,1)=theta0(i,1);
        x=zeros(length(v0),length(theta0),length(t));
        x(1)=0;
        y=zeros(length(v0),length(theta0),length(t));
        y(1)=2;
        for k=2:length(t)
            v(i,j,k)=v(i,j,k-1)+dt*(-g*sin(theta(i,j,k-1))-DL*(g/((vt)^2))*(v(i,j,k-1))^2);
            theta(i,j,k)=theta(i,j,k-1)+dt*(-(g/v(i,j,k-1))*cos(theta(i,j,k-1))+(g/((vt)^2))*v(i,j,k-1));
            x(i,j,k)=x(i,j,k-1)+dt*v(i,j,k-1)*cos(theta(i,j,k-1));
            y(i,j,k)=y(i,j,k-1)+dt*v(i,j,k-1)*sin(theta(i,j,k-1));
y(y(k)<0)=0;
        end
    end
end
surf(v0,theta0,x(:,:,end))

vmin=0;
vmax=15;
v0=linspace(vmin,vmax,100);
theta0=linspace(0,pi/2,100);
DL=1/5.2;
vt=5.5;
dt=0.1;
t=0:dt:10;
g=9.81;
xend=zeros(length(v0),length(theta0));
thetaend=zeros(length(v0),length(theta0));
vend=zeros(length(v0),length(theta0));
for i=1:length(v0)
    for j=1:length(theta0)
        v=zeros(1,length(t));
        v(1)=v0(i);
        theta=zeros(1,length(t));
        theta(1)=theta0(j);
        x=zeros(1,length(t));
        x(1)=0;
        y=zeros(1,length(t));
        y(1)=2;
        for k=2:length(t)
          
                v(k)=v(k-1)+dt*(-g*sin(theta(k-1))-DL*(g/((vt)^2))*(v(k-1))^2);
                theta(k)=theta(k-1)+dt*(-(g/v(k-1))*cos(theta(k-1))+(g/((vt)^2))*v(k-1));
                x(k)=x(k-1)+dt*v(k-1)*cos(theta(k-1));
                y(k)=y(k-1)+dt*v(k-1)*sin(theta(k-1));
                y(y(k)<0)=0;
        end
        vend(i,j)=v(end);
        thetaend(i,j)=theta(end);
        xend(i,j)=x(end);
    end
end
[v0,theta0]=meshgrid(linspace(vmin,vmax,100),linspace(0,2*pi,100));
surf(v0,theta0,xend)
%Plot the value of the final distance as a function of initial velocitty
%and angle.

%funtion called for part a
function dy=phugoid_ode(t,y)
DL=1/5.2;
vt=5.5;
g=9.81;
for i=1:length(t)
    dy(1,i)=-g*sin(y(2,i))-DL*(g/((vt)^2))*(y(1,i))^2;
    dy(2,i)=-(g/y(1,i))*cos(y(2,i))+(g/((vt)^2))*y(1,i);
    dy(3,i)=y(1,i)*cos(y(2,i));
    dy(4,i)=y(1,i)*sin(y(2,i));
end
end

The mnaxiimum value occured when v0=0 and theta naught equals two and forty three one hundredths radians.