Skip to content
No description, website, or topics provided.
MATLAB
Branch: master
Clone or download

Latest commit

Fetching latest commit…
Cannot retrieve the latest commit at this time.

Files

Permalink
Type Name Latest commit message Commit time
Failed to load latest commit information.
Grade
HeunsProb2.m
HeunsProb3.m
README.md
freefallheuns.m
freefallproblem.m
initvalprob2.m
initvalprob3.m
planeplot.png
prob2analytical.PNG
prob3analytical.PNG
prob4analytical.PNG
prob5me3255.m

README.md

06_initial_value_ode

#Problem 2 a) Analytically Solution 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 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

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

You can’t perform that action at this time.