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')
```