Newmark Method
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.
newmark(f, t0, t1, y0, ..., N = 100, zeta = 0.25, theta = 0.5)
f |
function in the differential equation y'' = f(x, y, y'); |
t0, t1 |
start and end points of the interval. |
y0 |
starting values as row or column vector;
|
N |
number of steps. |
zeta, theta |
two non-negative real numbers. |
... |
Additional parameters to be passed to the function. |
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.
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.
This is for demonstration purposes only; for real problems or applications
please use ode23
or rk4sys
.
Quarteroni, A., R. Sacco, and F. Saleri (2007). Numerical Mathematics. Second Edition, Springer-Verlag, Berlin Heidelberg.
# 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)
Please choose more modern alternatives, such as Google Chrome or Mozilla Firefox.