Skip to content
Permalink
e1c61d4854
Switch branches/tags

Name already in use

A tag already exists with the provided branch name. Many Git commands accept both tag and branch names, so creating this branch may cause unexpected behavior. Are you sure you want to create this branch?
Go to file
 
 
Cannot retrieve contributors at this time
200 lines (193 sloc) 3.8 KB
# 06_initial_value_ode
#Problem 2
a) Analytically
![Solution](./prob2analytical.PNG)
b) Eueler's Method
```
y_i=1;%initial conditions
n=3;%number of steps
h=1;%step size
y=zeros(1,n+1);
y(1)=y_i;
for i=1:n;
y(i+1)=y(i)-y(i)*h;
end
```
c)Heun's Approach
```
y_i=1;%initial conditions
n=3;%number of steps
h=1;%step size
y=zeros(1,n+1);
y(1)=y_i;
for i=1:n;
y_o=y(i)-y(i)*h;
y(i+1)=y(i)+((-y(i)+y_o)/2)*h;
end
```
| Time(s) | yTrue | yEuler | yHeuns |
| --- | --- | --- | --- |
| 0 | 1 | 1 | 1 |
| 1 | 0.3679 | 0 | 0.5 |
| 2 | 0.1353 | 0 | 0.25 |
| 3 | 0.0498 | 0 | 0.125 |
#Problem 3
a)
![Solution](./prob3analytical.PNG)
b)
```
y_i=1;%initial conditions
dy_i=0;
n=3;%number of steps
h=1;%step size
y=zeros(1,n+1);
y(1)=y_i;
for i=1:n;
y(i+1)=y(i)-3*(sin(3*(i-1)))*h;
end
```
c)
```
y_i=1;%initial conditions
n=3;%number of steps
h=1;%step size
y=zeros(1,n+1);
y(1)=y_i;
for i=1:n;
y(i+1)=y(i)+((-3*sin(3*(i-1))+-3*sin(3*(i))))/2*h;
end
```
| Time(s) | yTrue | yEuler | yHeuns |
| --- | --- | --- | --- |
| 0 | 1 | 1 | 1 |
| 1 | -0.99 | 1 | 0.7883 |
| 2 | 0.96 | 0.5766 | 0.9958 |
| 3 | -0.91 | 1.4149 | 0.7967 |
#Problem 4
a)
![Solution](./prob4analytical.PNG)
```
g=9.81;
cd=.25;
m=60;
for t=1:13
v(t)=sqrt((g*m)/cd)*tanh(sqrt((g*cd)/m)*(t-1));
x(t)=1000-(m/cd)*log(cosh(sqrt((g*cd)/m)*(t-1)));
end
```
b)
```
g=9.81;
cd=.25;
m=60;
x_i=1000;
v_i=0;
timesteps=12;
h=1;
v=zeros(1,timesteps+1);
x=zeros(1,timesteps+1);
v(1)=v_i;
x(1)=x_i;
for i=1:timesteps
v(i+1)=v(i)+(g-(cd/m)*v(i)^2)*h;
x(i+1)=x(i)-v(i)*h;
end
```
c)
```
g=9.81;
cd=.25;
m=60;
x_i=1000;
v_i=0;
timesteps=12;
h=1;
v=zeros(1,timesteps+1);
x=zeros(1,timesteps+1);
v(1)=v_i;
x(1)=x_i;
for i=1:timesteps
v_o=v(i)+(g-(cd/m)*v(i)^2)*h;
v(i+1)=v(i)+(((g-(cd/m)*v(i)^2)+v_o)/2)*h;
x_o=-v(i)*h;
x(i+1)=x(i)+(-v(i)+x_o)/2*h;
end
```
| Time (s)| yTrue (m)| yEuler (m) | yHeuns (m) |
| --- | --- | --- | --- |
| 0 | 1000 | 1000 | 1000 |
| 1 | 995.13 | 1000 | 1000 |
| 2 | 980.8924 | 990.19 | 990.190 |
| 3 | 958.323 | 970.97 | 966.066 |
| 4 | 928.8272 | 943.48 | 922.495 |
| 5 | 893.898 | 909.33 | 855.238 |
| 6 | 854.897 | 870.23 | 763.391 |
| 7 | 812.945 | 827.69 | 650.960 |
| 8 | 768.911 | 782.88 | 525.173 |
| 9 | 723.4515 | 736.62 | 392.609 |
| 10 | 671.962 | 689.47 | 257.175 |
| 11 | 629.82 | 641.78 | 120.640 |
| 12 | 582.22 | 593.75 | -16.298 |
| Time Step | vTrue | vEuler | vHuens |
| --- | --- | --- | --- |
| 0 | 0 | 0 | 0 |
| 1 | 9.68 | 9.81 | 9.81 |
| 2 | 18.61 | 19.21902 | 24.12402 |
| 3 | 26.28 | 27.48997 | 43.57116 |
| 4 | 32.45 | 34.15123 | 67.25655 |
| 5 | 37.17 | 39.10162 | 91.84714 |
| 6 | 40.64 | 42.54105 | 112.43114 |
| 7 | 43.11 | 44.81046 | 125.78687 |
| 8 | 44.84 | 46.25389 | 132.56390 |
| 9 | 46.07 | 47.14963 | 135.43424 |
| 10 | 46.84 | 47.69676 | 136.53455 |
| 11 | 47.59 | 48.02768 | 136.93815 |
| 12 | 47.77 | 48.22660 | 137.08365 |
#Problem 5
```
g=9.8;
v_t=5.5;
Cd=1/40;
Cl=1;
v_i=10;
theta_i=0;
y_i=2;
x_i=0;
t=20;
h1=0.1;
h2=0.01;
dt1=t/h1;
dt2=t/h2;
v1=zeros(1,dt1+1);
v1(1)=v_i;
theta1=zeros(1,dt1+1);
theta1(1)=theta_i;
y1=zeros(1,dt1+1);
y1(1)=y_i;
x1=zeros(1,dt1+1);
x1(1)=x_i;
v2=zeros(1,dt2+1);
v2(1)=v_i;
theta2=zeros(1,dt2+1);
theta2(1)=theta_i;
y2=zeros(1,dt2+1);
y2(1)=y_i;
x2=zeros(1,dt2+1);
x2(1)=x_i;
for i=1:dt1
v1(i+1)=v1(i)+h1*(-g*sin(theta1(i))-(Cd/Cl)*(g/v_t^2)*(v1(i)^2));
theta1(i+1)=theta1(i)+h1*(-(g/v1(i))*cos(theta1(i))+(g/v_t^2)*v1(i));
x1(i+1)=x1(i)+h1*v1(i)*cos(theta1(i));
y1(i+1)=y1(i)+h1*v1(i)*sin(theta1(i));
end
for i=1:dt2
v2(i+1)=v2(i)+h2*(-g*sin(theta2(i))-(Cd/Cl)*(g/v_t^2)*(v2(i)^2));
theta2(i+1)=theta2(i)+h2*(-(g/v2(i))*cos(theta2(i))+(g/v_t^2)*v2(i));
x2(i+1)=x2(i)+h2*v2(i)*cos(theta2(i));
y2(i+1)=y2(i)+h2*v2(i)*sin(theta2(i));
end
plot(x1,y1,'l',x2,y2)
xlabel('x (m)')
ylabel('y (m)')
```
![Plane Flight](./planeplot.png)