Stetson University Create a Python Code Computer Science Task

Assignment #6BME 3420 – Fall 2021
Dec 02, 2021
1- [5 points]
Another approach to solve an ODE using the Explicit Euler formula, is to use the formula to directly
find an equation that predicts S(t
j+1
)
given S(t ). For example, for the linearized pendulum
j
problem, where
0
dS(t)
= [
dt
−
1
] S(t)
g
0
l
we can use the Explicit Euler Formula to obtain
0
S(tj+1 ) = S(tj ) + h [
1
S(tj+1 ) = [
0
1
g
−
0
l
0
0
1
] S(tj ) + h [
1
S(tj+1 ) = [
1
g
−
l
0
] S(tj )
h
] S(tj )
gh
−
] S(tj )
1
l
However, this solution has the same problem as before. Here is a code that implements this solution
for the pendulum problem with with
1
S(0) = [
]
0
and
g
√
= 4
l
In [6]:
import numpy as np
import matplotlib.pyplot as plt
plt.style.use(‘seaborn-poster’)
h=0.1
time = np.arange(0,5+h,h)
s0 = np.array([[1],[0]])
s = np.zeros((len(time),2))
s[0,:] = s0.T
w = 4 #(g/l)^(1/2)
for i in range(1,len(time)-1):
s[i,:] = np.array([[1, h],[-(w**2)*h, 1]])@s[i-1,:]
plt.figure(figsize = (5, 5))
plt.plot(time, np.cos(w*time),label=’Exact Solution’)
plt.plot(time,s[:,0], label = ‘Explicit Euler Formula’)
plt.xlabel(‘Time (s)’)
plt.ylabel(‘$\Theta(t)$’)
plt.legend()
plt.show()
An alternative is to implement the Intrinsic Euler Formula, which approximate the state S(t
dS(tj+1 )
S(tj+1 ) = S(tj ) + (tj+1 − tj )
Note that this equation depends on
dS(tj+1 )
dt
,
dt
which is unknonw!
However, we can use some algebra to solve this problem and obtain
dS(tj+1 )
0
= [
dt
1
g
−
Reemplazing in the Implicit Euler formula we obtain
l
0
] S(tj+1 )
j+1 )
by
0
S(tj+1 ) = S(tj ) + h [
[
1
0
0
−
1
[
g
l
0
l
] S(tj+1 ) = S(tj )
−h
h
1
] S(tj+1 ) = S(tj )
S(tj+1 ) = [
for any S(t
]
gh
j
−1
−h
S(tj )
1
l
)
] S(tj+1 ),
1
1
j+1
0
l
g
1
Which allow us to compute S(t
−
0
] S(tj+1 ) − h [
1
g
)
a) [2 points] Implement an algorithm that allows you to solve the pendulum problem using the
Intrisic Euler Formula. Plot the numerical and exact predictions, discuss your plot.
Hint : take a look at the code for the Extrinsic Euler Formula and make the necesary changes
As you observed in point a, the Intrisic Euler Formula does not provide a correct solution either. One
approach overshoots the solution and another undershoots the solution. A reasonable approach
then, is to find the average between the Intrinsic and Extrinsic formulas, this is called the trapezoidal
formula, and is given by
dS(tj )
h
S(tj+1 ) = S(tj ) +
[
dS(tj+1 )
+
dt
2
],
dt
b) [1 point] Demonstrate that for the pendulum problem, where
dS(tj )
0
= [
dt
−
dS(tj+1 )
the trapezoidal formula is given by
l
0
= [
dt
1
g
0
1
g
−
] S(tj )
l
0
] S(tj+1 )
⎡ 1
S(tj+1 ) =
⎣
−
gh
2
⎤
⎦
1
2l
−1
h
⎡
h
1
⎣−
2
gh
⎤
S(tj )
1 ⎦
2l
You should write the mathematical operations to derive this formula using pen and paper.
c) [2 points] Implement an algorithm that allows you to solve the pendulum problem using the
Trapezoidal Formula. Plot the numerical and exact predictions, discuss your plot.
2- [5 points]
A model that describes the relation between position, velocity, and acceleration of a mass in a
simple sprin-mass-damper system is given by
mẍ + bẋ + kx = F (t),
where F (t) is the force applied to the mass.
a) [3 points] Assume that F (t)
= 0
when b
. Create a plot for each condition and discuss your observations.
,
= 0 b = 1
, and b
= 10
For this problem, assume that m
for all t. Solve the ODE to determine the position of the mass
= 1
and k
= 5
. The initial conditions are
x(0) = 1
ẋ(0) = 0
and you are interested in the solution between [0, 20] with h
= 0.01
Hint: Remember that numerical methods can only deal with first order equations, so you need to
reduce the ODE order before solving the problem.
b) [2 points]
⎧0
F (t) = ⎨ 1
⎩
0
if t < 3 if 4 ≤ t ≤ 8 otherwise Solve the ODE to determine the position of the mass when b condition and discuss your observations. = 1 and b = 10 . Create a plot for each For this problem, assume that m = 1 and k = 5 . The initial conditions are x(0) = 1 ẋ(0) = 0 and you are interested in the solution between [0, 20] with h In [ ]: = 0.01