newmark: Newmark Method

Description

Newmark's is a method to solve higher-order differential equations without passing through the equivalent first-order system. It generalizes the so-called `leap-frog' method. Here it is restricted to second-order equations.

Usage

newmark(f, t0, t1, y0, ..., N = 100, zeta = 0.25, theta = 0.5)

Arguments

function in the differential equation $y'' = f(x, y, y')$; defined as a function $R \times R^2 \rightarrow R$.

t0, t1

start and end points of the interval.

starting values as row or column vector; y0 needs to be a vector of length 2, the first component representing y(t0), the second dy/dt(t0).

number of steps.

zeta, theta

two non-negative real numbers.

...

Additional parameters to be passed to the function.

Value

List with components t for grid (or `time') points between t0 and t1, and y an n-by-2 matrix with solution variables in columns, i.e. each row contains one time stamp.

Details

Solves second order differential equations using the Newmark method on an equispaced grid of N steps.

Function f must return a vector, whose elements hold the evaluation of f(t,y), of the same dimension as y0. Each row in the solution array Y corresponds to a time returned in t.

The method is `implicit' unless zeta=theta=0, second order if theta=1/2 and first order accurate if theta!=1/2. theta>=1/2 ensures stability. The condition set theta=1/2; zeta=1/4 (the defaults) is a popular approach that is unconditionally stable, but introduces oscillatory spurious solutions on long time intervals. (For these simulations it is preferable to use theta>1/2 and zeta>(theta+1/2)^(1/2).)

No attempt is made to catch any errors in the root finding functions.

References

Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.

Examples

Run this code

# Mathematical pendulum  m l y'' + m g sin(y) = 0
pendel <- function(t, y)  -sin(y[1])
sol <- newmark(pendel, 0, 4*pi, c(pi/4, 0))

## Not run: 
# plot(sol$t, sol$y[, 1], type="l", col="blue",
#      xlab="Time", ylab="Elongation/Speed", main="Mathematical Pendulum")
# lines(sol$t, sol$y[, 2], col="darkgreen")
# grid()## End(Not run)

Run the code above in your browser using DataCamp Workspace