program留学生代做、Python编程设计调试、代写data、代写Python

日期：2020-05-02 11:23

Page 1 of 9

Assignment 9: Explicit Time-Marching Schemes

Due Date: Wednesday 29 April 11:59pm.

Overview

This assignment brings together object-oriented programing with numerical time-marching

methods. For the first part of the assignment you will implement some explicit time-marching

schemes within an object-oriented framework to solve systems of ODEs that model initial value

problems. For the second part of the assignment you will model the dynamics of a crane cable

that might fail when the crane mechanism seizes up while lowering a load.

Instructions and Deliverables

Part 1: ode module

Create a module file named ode.py. This file must define the following classes:

- ODESolver: a base class for all time-marching methods.

- ForwardEuler: a subclass of ODESolver that implements forward Euler.

- RK4: a subclass of ODESolver that implements the 4th-order Runge-Kutta method.

Before describing the methods to implement for each of the above classes, we first illustrate how

the classes would be used to solve an initial-value problem (IVP). Suppose we want to solve the

bacteria concentration IVP described in the lectures:

This could be solved and the solution displayed using the ForwardEuler class with the

following code:

import ode

import numpy as np

def bacteria(c, t):

"""Returns the ODE right-hand-side function

for the bacteria decontamination model"""

f = -0.6*c

return f

tspan = np.array([0.0, 7.0])

num_steps = 10

y0 = np.array([1.0e7]) # note that the IC is passed as an array

solver = ode.ForwardEuler(bacteria) # construct variable of type ForwardEuler

t, y = solver.solve_ivp(tspan, y0, num_steps) # uses method solve_ivp

Page 2 of 9

print()

print(" i time count ")

print("--- ---- --------")

for i in range(y.shape[0]):

string = '{:3d}'.format(i ) + " "

string += '{:4.2f}'.format(t[i] ) + " "

string += '{:8.2e}'.format(y[0, i])

print(string)

Look closely at the follwing two lines in particular. You must make yourself familiar with this

object oriented approach to applying the numerical methods to understand this assignment.

solver = ode.ForwardEuler(bacteria) # construct variable of type ForwardEuler

– ode is the imported module file that contains various numerical methods for solving

ODEs.

– solver is a variable constructed to be of type ForwardEuler, a class that is

defined in module ode. solver gets its properties from class ForwardEuler

(and its super class ODESolver).

– bacteria is a function (defined above) attached to solver via the __init__

function that ForwardEuler will inherit from super class ODESolver.

Through this construction, variable solver now carries bacteria with it to use in the

solution method.

t, c = solver.solve_ivp(tspan, c0, num_steps) # uses method solve_ivp

– solver is the variable constructed above with the function bacteria attached.

– solve_ivp is a method in class ForwardEuler (that it actually inherits from

super class ODESolver) which will apply the forward Euler numerical method

to solve the initial value problem. solve_ivp takes three arguments:

0

B tspan is a two element NumPy array holding tspan[0] = t and

tspan[1] f

= t .

B y0 is a 1D NumPy array with as many elements as there are ODEs in

the system of ODEs to be solved. As shown above, the bacteria

example has only one ODE so a one-element array holds the one

initial condtion.

0 f B num_steps is the number of time steps to use between t and t .

To solve the same IVP with the 4th-order Runge-Kutta method, we would only need to replace

the ForwardEuler constructor with the RK4 constructor.

Page 3 of 9

The ODESolver base class shall meet the following requirements:

C The class constructor will take in a function and define the _func member based on

this function. Thus, the first-line definition should be

def __init__(self, func):

For any variable constructed by ODESolver we want the argument passed to

__init__ to be “attached” to that constructed variable.

Recall from Lecture 16 how __init__ works. In the example on pp. 6-11,

when we construct a variable of type Vector we associate parameters in

function __init__ (x and y) by which arguments can be assigned to the

constructed variable (self._x and self._y). In ODESolver we do the

same, except that we associate a parameter that can be a function: func, in the

same way that we can pass functions to other functions (e.g., see Lecture 19 p. 14

where the forward_euler function has parameter func). Note that we want

to use an underscore as the first character to signal that the attached function is

private: that when interacting with your class no one should arbitrarily reassign

self._func for any reason. (The same way one should not write something

like math.pi = 3; perhaps it should instead be named math._pi to signal

that.)

C The class should include the member function solve_ivp with the following

first-line definition:

def solve_ivp(self, tspan, y0, num_steps):

This was described in detail above.

As with the in-class exercises, raise a ValueError with a useful message if the

arguments to tspan or num_steps contain an obvious error.

The method solve_ivp should return two arrays, t and y, in that order. The

array t should hold all num_steps+1 time values including tspan[0] and

tspan[1]. The array y should be a two-dimensional NumPy array with

y.shape[0] = y0.size and y.shape[1] = t.size; that is, the

number of “rows” in y should equal the number of dependent variables (and, thus,

the number of ODEs in the system), and the number of “columns” in y should

equal the number of values in t.

Finally, solve_ivp should not be re-implemented in the subclasses. The idea

is to write this routine once and use it for all the subclasses via inheritance. The

Page 4 of 9

solution will thus be calculated by solve_ivp as

for n in range(0, num_steps):

y[:,n+1] = self.step(t[n], dt, y[:,n])

Compare this to how solutions were calculated in previous in-class exercises (e.g.,

using forward Euler, and note two changes:

– the solution variable, y, is now a 2D NumPy array

– the step method comes from self. For the bacterial colony example

above, self is the forwardEuler type that solver was

constructed to be. Thus, self.step will be the same as calling

forwardEuler.step – the step function defined in the

forwardEuler sub-class. (If the user instead chose to construct

solver to be of type RK4, then this would be the same as calling

RK4.step and would invoke the step function defined in that

sub-class.)

The ODESolver class is an abstract base class: it should be used to derive new

classes such as ForwardEuler and RK4, and it should not be used to solve

IVPs on its own. To enforce this intent, provide ODESolver with a step method

with the following first-line definition:

def step(self, t_old, dt, y_old):

However, for the ODESolver base class only, this function should only raise a

NotImplementedError and not do anything else.

The ForwardEuler subclass should meet the following requirements.

C The ForwardEuler class should be derived from ODESolver.

C The ForwardEuler class should implement only one method, namely step.

The step method should have the identical first-line definition as for super class

ODESolver:

def step(self, t_old, dt, y_old):

where t_old n

is the time value at the current time (referred to as t in the lecture

slides), dt is the time step size )t, and y_old is the variable value at the current

n n+1 time (y ). step will return the value y as defined by the forward Euler

numerical method.

Page 5 of 9

The RK4 subclass should meet the following requirements.

C The RK4 class should be derived from ODESolver.

C The RK4 class should implement only one method, namely step. The step method

should have the first-line definition shown above for forwardEuler, and it

n+1 should return the value y as defined by the classical 4th-order Runge-Kutta

method.

As usual, provide docstrings for the ode module, the classes, and all methods.

Part 2: cable mechanism seizure failure

Background and Models

Consider the classic crane from IEA, pictured at right. A force

balance on the load when static gives:

where

+y is in the downward (lowering) direction

W is the weight of the load [N]

T is the tension in the cable [N]

m is the mass of the load [kg] and g is the constant of gravitation [9.80665 m/s ]

2

Recall that “static” means non-accelerating, so when a load is being lowered at a constant

velocity it is also in static equilibrium. From this model, a cable that can support a 7 kN tension

can theoretically lower a 7 kN load at constant velocity.

However, what if the cable mechanism suddenly seizes while lowering the load? Suppose that

the load is being lowered at some velocity v and suddenly the cable mechanism stops. From

Physics I, the dynamics of the load involves:

– inertia (the “ma” in F = ma),

– a damping/restoring force (not present in this problem as described so far) proportional

to the dynamic velocity of the load,

– the elastic tension in the cable (a spring force) due to elongation from the change in

momentum of the load when the cable seizes,

– the static tension in the cable, and

– the weight of the load,

as shown in this equation:

Page 6 of 9

where

a(t) is the acceleration of the load [m/s ]

2

v(t) is the velocity of the load [m/s]

y(t) is the displacement of the load from the point it stops [m]

b is a damping coefficient [N@s/m]

k is the spring constant of the elastic cable (stretching from the inertia) [N/m]

The initial conditions are that the initial position y(0) is arbitrarily set to zero, and the initial

velocity v(0) is the ve initial locity at which the load is being lowered, v .

From the static model W ! T = 0. Noting that dy(t)/dt = v(t) and d y(t)/dt = a(t), and 2 2

multiplying through by !1, the equation of motion of the dynamic load can be written as:

After some algebra this can also be written as a system of first order ODEs:

In Physics I and possibly in Introduction to Differential Equations you have already seen these

derived analytically, which is useful for checking a numerical solution. For a review of the

theory of harmonic oscillators you may find the following helpful:

https://physics.info/sho/

http://hyperphysics.phy-astr.gsu.edu/hbase/shm.html#c1

https://en.wikipedia.org/wiki/Harmonic_oscillator

However, you may not use the analytical solutions for this assignment (except as optional

checks). All results should be computed from your numerical solution.

Frequency f of undamped (b = 0) harmonic oscillator

Note: if you plot y(t) vs t for 1 second you should be able to count f periods if

you’ve solved the system correctly.

Page 7 of 9

Max max imum tension T in the cable

Max max imum stress F in the cable

Criterion for cable failure

Deliverables

In a module file titled hw9.py include the following (as well as the usual useful doctstrings and

comments throughout).

final C all model input data (e.g., mass, spring constant, etc.) and numerics input (t , number

of steps) should be included as module variables near the top of your file that you

can easily adjust. Include the units for the variables as in-line comments. For

example,

g = 9.80665 # constant of gravitation [m/s**2]

You may find it helpful to optionally print this out for your information; e.g.,

print()

print("LOAD ----- ----- -----")

print(" mass, m [kg] = ", '{:7.1f}'.format(m))

print(" weight, W [N] = ", '{:7.1f}'.format(W))

print("CABLE ---- ----- -----")

print(" spring constant, k [N/m] = ", '{:12.2e}'.format(k))

# etc.

C a function that defines the right hand side of the ODE model to solve.

This function will input the array of state variables and the time and output the

right hand side of the system of ODEs model; e.g.,

Page 8 of 9

def crane_cable(state, t):

'''Returns the ODE right-hand-side function

for the crane cable mechanism seizure failure model

Input: state[:,i] - 1D slice of state variables at i-th time step

t - time at i-th time step

Returns: dsdt[:,i] - 1D ODE rhs function evaluated at input state

'''

C code to import and apply an appropriate ODE solver from your ode.py module.

C code to create a professional .png plot of load position versus time for 1 second. As

described above, use this plot to verify that you have solved the model correctly.

Use the if(Make_Plot): technique described in Assignment 8 to hide any

references to matplotlib from Submitty before you submit file hw9.py.

C code to compute the principle “underlined” quantities listed above and print them out

in a well formatted table that includes units. When you are ready to submit your

assignment to Submitty, copy this output into comments at the end of file hw9.py.

There will be two sets of output, described subsequently. The following is an

example of how you might do this. Please output the values in the units shown.

print("RESULTS ** ***** ***** ***** *****")

print(" frequency =", '{:6.2f}'.format(frequency) , "[Hz]")

print(" amplitude =", '{:6.2f}'.format(max_amplitude*1.0e3), "[mm]")

print(" tension =", '{:6.2f}'.format(max_tension/1.0e3) , "[kN]")

print(" stress =", '{:6.2f}'.format(max_stress /1.0e6) , "[MPa]")

C a comment block (e.g., text with000 A)... B)... C)... D)...000

at the very end of your code briefly answering the following questions.

A) How do you know you’ve used an appropriate ODE solver?

[Hint: look over your notes that describe the stability of the methods you

have available in your ode.py module.]

B) State a criterion for deciding how many time steps to use (i.e., how to get a

small enough )t) and state the number of time steps you used to meet that

criterion. You do not necessarily need a quantitative criterion.

[Hint: our derivations in the notes were about the stability of the solution,

not the accuracy. How can you judge accuracy? Consider the

convergence criteria we’ve studied previously this semester, and run

Page 9 of 9

numerical experiments with various numbers of steps to come to a

conclusion.]

C) Choose a unique cable strength between 50.0 and 100 MPa. Use your code to

manually (or automatically, if you like) compute the diameter of the cable

(in mm, to two significant figures rounded up) needed to survive this

accident with a factor of safety of 3.0. (Note that factors of safety are still

commonly used in design scoping work, but are now deprecated in many

safety related engineering analyses.)

[Hint: “manually” means you can manually adjust an input value and

rerun your code to get a desired output. Ideally you would implement a

root-finding procedure to do this, but unfortunately (?) we believe the time

remaining in this semester is too short to ask for that.]

D) A vendor has a damping bar attachment that can mitigate a cable seizure

accident with a damping coefficient of b = 8000.0 N@s/m. For the cable

strength you chose in Question C, qualitatively how does this change the

cable diameter – a lot or a little, and why?

[Hint: it is recommended that you carefully examine and compare the two

plots and recall how the maximum stress is calculated to consider this

question.]

Please note that you have a lot of freedom in how you code and answer this assignment. The

point is not to get a “right” answer, it is to identify a safe design that you are confident it without

knowing exactly what the “right” answer is.

Use the following data in your model. Recall that these should apper as module variables.

mass of load [kg] = 500.0 kg

damping coefficient [N@s/m] = 0.0 N.s/m (except for question D)

spring constant of cable in N/m = 1.05@10 N/m

6

constant of gravitation = 9.80665 m/s

cable speed = 0.1 m/s (This assignment should be

informative about why the speed is

generally so slow; experiment with

other values to see the effect.)

final t = 1.0 s