Homework answers / question archive / Assignment #6 BME 3420 - Fall 2021 Another approach to solve an ODE using the Explicit Euler formula, is to use the formula to directly find an equation that predicts S(tj+1) given S(tj)

Assignment #6 BME 3420 - Fall 2021 Another approach to solve an ODE using the Explicit Euler formula, is to use the formula to directly find an equation that predicts S(tj+1) given S(tj)

Computer Science

Share With

---

Assignment #6

BME 3420 - Fall 2021

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_j). For example, for the linearized pendulum problem, where

dt = [ − g 0]S(t) dS(t) 0 1l

we can use the Explicit Euler Formula to obtain

1. 0. [

      1

                                                                                   S(t_j+1) = S(t_j) + h_g                ]S(t_j)

                                                                                                                     − _l       0

1. 1. 0  0  1

                                                                    S(t_j+1) = [                  ]S(t_j) + h[         g                                                                                          ]S(t_j)

0. 1    − _l         0

1. [

   h

                                                                                         S(tj+1) =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

√

_l

 

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()

In [6]:

 

 

An alternative is to implement the Intrinsic Euler Formula, which approximate the state S(t_j+1) by

dS(t_j+1)

S(t_j+1) = S(t_j) + (t_j+1 − t_j) , dt

dS(t_j+1)

Note that this equation depends on  which is unknonw! dt

However, we can use some algebra to solve this problem and obtain

dS(tj+1)                 0  = [         g        dt                  − l

1 S(t 0]           j+1)

Reemplazing in the Implicit Euler formula we obtain

 

[

                                                                                                                     0       1

 

                                                                                S(t_j+1) = S(t_j) + h_g              ]S(t_j+1),

− 0 l

                                                                  1     0                                  0       1

                                                                                 ]S(t_j+1) − h[          g         ]S(t_j+1) = S(t_j)

0. 1 − _l 0

1. [

   −h h 1

                                                                                      g                          ]S(t_j+1) = S(t_j)

l

 

[

                                                                                                          1       −h     −1

 

                                                                                      S(tj+1) =gh                        ]    S(tj)

1 l

Which allow us to compute S(t_j+1) for any S(t_j)

1. [2 points] Implement an algorithm that allows you to solve the pendulum problem using theIntrisic 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

                                                                                                       h      dS(t_j) dS(t_j+1)

                                                                      S(t_j+1) = S(t_j) +                                                            [                 + ],

                                                                                                       2         dt            dt

2. [1 point] Demonstrate that for the pendulum problem, where

dS(tj)               0        = [       g          dt                − l	1 ]S(t ) 0           j
dS(tj+1)                 0  = [         g        dt                  − l	1 S(t 0]           j+1)

the trapezoidal formula is given by

 

?

1

h

2

g

h

2

l

S(tj+1) =

?

 gh          2 ??        ??        ??S(tj) 1        − h          −1          1

 

                                                                                                  1                 −

₂l

You should write the mathematical operations to derive this formula using pen and paper.

3. [2 points] Implement an algorithm that allows you to solve the pendulum problem using theTrapezoidal 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

mx¨ + bx? + kx = F(t),

where F(t) is the force applied to the mass.

1. [3 points] Assume that F(t) = 0 for all t. Solve the ODE to determine the position of the mass when b = 0, b = 1, and b = 10. Create a plot for each condition and discuss your observations.

For this problem, assume that m = 1 and k = 5. The initial conditions are

x(0) = 1 x?(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.

2. [2 points]

                                                                                                 ???0    if t < 3

                                                                                      F(t) =        1       if 4 ≤ t ≤ 8

                                                                                                        0     otherwise

Solve the ODE to determine the position of the mass when b = 1 and b = 10. Create a plot for each condition and discuss your observations.

For this problem, assume that m = 1 and k = 5. The initial conditions are

x(0) = 1 x?(0) = 0

and you are interested in the solution between [0,20] with h = 0.01

In [ ]: