ODE Solu)on Methods
§ Euler Method (1st order) x n +1 = x n + f (x n )Δt
§ Improved Euler Method (2nd order) x˜ n +1 = x n + f (x n )Δt
1
x n +1 = x n + [ f (x n ) + f ( x˜ n +1 )]Δt
2
§ Runge-‐Ku>a Methods (higher order) k1 = f (x n )Δt,
x n +1 = x n +
k 2 = f (x n +
1
k1 )Δt,
2
1
[k1 + 2k 2 + 2k3 + k 4 ]Δt
6
k 3 = f (x n +
1
k2 )Δt,
2
k 4 = f (x n + k 3 )Δt
Newton’s Law Problems
§ Second order ODE
§ Constant Force F=C (kinema)cs)
§ Time Dependent Force F=F(t)
§ Velocity Dependent Force F=F(v) (ex. Projec)le with drag)
§ Posi)on Dependent Force F=F(x) (ex. Simple harmonic mo)on)
§How can ODE solu)on methods help?
€
Newton’s Law Problems
§ Second order ODE
∑ F = ma
§ Constant Force F=C (kinema)cs)
§ Time Dependent Force F=F(t)
§ Velocity Dependent Force F=F(v) (ex. Projec)le with drag)
§ Posi)on Dependent Force F=F(x) (ex. Simple harmonic mo)on)
§How can ODE solu)on methods help?
€
Newton’s Law Problems
§ Second order ODE
∑ F = ma
§ Constant Force F=C (kinema)cs)
§ Time Dependent Force F=F(t)
§ Velocity Dependent Force F=F(v) (ex. Projec)le with drag)
§ Posi)on Dependent Force F=F(x) (ex. Simple harmonic mo)on)
§How can ODE solu)on methods help?
F = −kx = ma
˙x˙ = −ω 2 x
€
Newton’s Law Problems
§ Second order ODE
∑ F = ma
§ Constant Force F=C (kinema)cs)
§ Time Dependent Force F=F(t)
§ Velocity Dependent Force F=F(v) (ex. Projec)le with drag)
§ Posi)on Dependent Force F=F(x) (ex. Simple harmonic mo)on)
§How can ODE solu)on methods help?
F = −kx = ma
˙x˙ = −ω 2 x
x˙ = v
v˙ = −ω 2 x
Damped SHM
clear all
%Shows how to solve the damped SHM problem
for three %values of the damping constant
C(1)=1;
C(2)=50;
C(3)=200;
for dampcons_choice=1:3
m=20; % m=mass
c=C(dampcons_choice); %c=damping
factor
k=20; %spring constant
final_time=100;
initial_time=0;
delta_t=.01;
N=final_time/delta_t;
v=zeros(N,1);
x=zeros(N,1);
t=zeros(N,1);
v(1)=0;
x(1)=1;
t(1)=initial_time;
Damped SHM
for i=1:N-1
xk_1=v(i)*delta_t;
vk_1=((-c/m)*v(i)-(k/m)*x(i))*delta_t;
xk_2=(v(i)+.5*vk_1)*delta_t;
vk_2=((-c/m)*(v(i)+.5*vk_1)-(k/m)*(x(i)
+xk_1))*delta_t;
xk_3=(v(i)+.5*vk_2)*delta_t;
vk_3=((-c/m)*(v(i)+.5*vk_2)-(k/m)*(x(i)
+xk_2))*delta_t;
xk_4=(v(i)+.5*vk_3)*delta_t;
vk_4=((-c/m)*(v(i)+.5*vk_3)-(k/m)*(x(i)
+xk_3))*delta_t;
x(i+1)=x(i)
+(1/6)*(xk_1+2*xk_2+2*xk_3+xk_4);
v(i+1)=v(i)
+(1/6)*(vk_1+2*vk_2+2*vk_3+vk_4);
t(i+1)=t(i)+delta_t;
end
Damped SHM
figure(1)
subplot(2,1,1)
hold on
if(dampcons_choice==1)
plot(t,x,'b')
elseif(dampcons_choice==2)
plot(t,x,'r')
elseif(dampcons_choice==3)
plot(t,x,'k')
end
subplot(2,1,2)
hold on
if(dampcons_choice==1)
plot(t,v,'b')
elseif(dampcons_choice==2)
plot(t,v,'r')
elseif(dampcons_choice==3)
plot(t,v,'k')
end
end
ODE Solu)on Methods
§ Adap)ve Step Size
§ step halving
§ embedded RK methods
§ ODE S)ffness
§ implicit vs explicit solns
dy
= −ay
dt
dy
= −ay
dt
dy i+1
h
dt
y i+1 = y i − ay i+1h
yi
y i+1 =
1+ ah
dy
y i+1 = y i + i h
dt
y i+1 = y i − ay i h
y i+1 = y i (1− ah)
€
y i+1 = y i +
€
€
ODE Solu)on Methods
§ Adap)ve Step Size
§ step halving
§ embedded RK methods
§ ODE S)ffness
§ implicit vs explicit solns
dy
= −1000y + 3000 − 2000e−t
dt
y = 3 − 0.998e−1000t − 2.002e−t
ODE Solvers in MATLAB
>> help ode45
ODE45 Solve non-‐s)ff differen)al equa)ons, medium order method.
[TOUT,YOUT] = ODE45(ODEFUN,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system of differen)al equa)ons y' = f(t,y) from )me T0 to TFINAL with ini)al condi)ons Y0. ODEFUN is a func)on handle. For a scalar T
and a vector Y, ODEFUN(T,Y) must return a column vector corresponding to f(t,y). Each row in the solu)on array YOUT corresponds to a )me returned in the column vector TOUT. To obtain solu)ons at specific )mes T0,T1,...,TFINAL (all increasing or all decreasing), use TSPAN = [T0 T1 ... TFINAL].
ODE Solvers -‐ Bungee with Cord
§ Solve the bungee jumper problem including effect of cord
§ system of 2 ode’s
dx
=v
dt
€
ODE Solvers -‐ Bungee with Cord
§ Solve the bungee jumper problem including effect of cord
§ system of 2 ode’s
dx
=v
dt
§ and
§if cord is slack
€
dv
c
= g − sign(v) d v 2
dt
m
€
ODE Solvers -‐ Bungee with Cord
§ Solve the bungee jumper problem including effect of cord
§ system of 2 ode’s
dx
=v
dt
§ and
§if cord is slack
§if cord is taught
€
dv
c
= g − sign(v) d v 2
dt
m
€
dv
c
k
γ
= g − sign(v) d v 2 − (x − L) − v
dt
m
m
m
€
ODE Solvers -‐ Bungee with Cord
§ Solve the bungee jumper problem including effect of cord
§ system of 2 ode’s
dx
=v
dt
§ and
§if cord is slack
§if cord is taught
€
dv
c
= g − sign(v) d v 2
dt
m
€
dv
c
k
γ
= g − sign(v) d v 2 − (x − L) − v
dt
m
m
m
€
ODE Solvers -‐ Bungee with Cord
§ Solve the bungee jumper problem including effect of cord
§ system of 2 ode’s
dx
=v
dt
§ and
§if cord is slack
§if cord is taught
€
dv
c
= g − sign(v) d v 2
dt
m
€
dv
c
k
γ
= g − sign(v) d v 2 − (x − L) − v
dt
m
m
m
€
ODE Solvers -‐ Bungee with Cord
§ Solve the bungee jumper problem including effect of cord
function dydt = bungeecord(t,y,L,cd,m,k,gamma)
g=9.81;
cord=0;
if(y(1)>L) %determine if the cord exerts a force
cord=k/m*(y(1)-L)+gamma/m*y(2);
end
dydt=[y(2); g-sign(y(2))*cd/m*y(2)^2-cord];
[t,y]=ode45(@bungeecord,[0 50],[0 0],[],30,0.25,68.1,40,8);
plot(t,-‐y(:,1),'-‐',t,y(:,2),':')
legend('x(m)','v(m/s)')
Spinning Tennis Ball Problem
z
v
ω
D
G
M
velocity
ω
angular velocity
drag
gravity
€
magnus
v
D
G
€
€
x
€
€
M
y
Spinning Tennis Ball Problem
Spinning Tennis Ball Problem
HW #2 – Use MATLAB to solve for the posi)on of a tennis ball through )e from the moment
it is given an ini)al velocity of Vo un)l it hits the ground for three cases
case 1) in a vacuum
case 2) in air but not spinning
case 3) in air and spinning with angular velocity ω
As a final product of this solu)on plot the trajectory for the three cases.
h=1m, ω=20, Vo=25 m/s, θ=15 degrees, g=9.8 m/s^2, d=.063m, m=.05kg, rho 1.29
© Copyright 2024 ExpyDoc
