-
-
Save anonymous/1442252 to your computer and use it in GitHub Desktop.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
## a second order differential equation | |
## D2 y + 2 r w D y + w^2 y = 0 | |
# D is the differential operator | |
# r is the damping ratio (r = c/(2mw)) | |
# w is the undamped frequency | |
# c is the viscous damping coefficent | |
# m is the mass | |
#//////////////////////////////////////////////////////////////////////////////// | |
#// // | |
#// Description: // | |
#// A differential equation of order greater than one can be transformed // | |
#// to a multidimensional first order differential equation. In // | |
#// particular, given a second order differential equation y'' = f(x,y,y') // | |
#// define z = y', then the pair (y,z)' = (z, f(x,y,z) ) is a first order // | |
#// differential equation. The initial conditions y(x0) = y[0], // | |
#// y'(x0) = c becomes the initial condition (y,z)(x0) = (y[0],c). // | |
#// when x = x0 numerically evaluates f(x,y) four times per step. Then for // | |
#// step i+1, // | |
#// y[i+1] = y[i] + h * ( yp[i] + 1/6 (k1 + k2 + k3 ) and // | |
#// yp[i+1] = yp[i] + 1/6 * ( k1 + 2 k2 + 2 k3 + k4 ) where // | |
#// k1 = h * f(x[i], y[i], yp[i]), // | |
#// k2 = h * f(x[i]+h/2, y[i]+(h/2)*yp[i], yp[i]+k1/2), // | |
#// k3 = h * f(x[i]+h/2, y[i]+(h/2)*yp[i]+(h/4)*k1, yp[i]+k2/2), // | |
#// x[i] = x0 + i * h. // | |
#// // | |
#//////////////////////////////////////////////////////////////////////////////// | |
f <- function(t_n,y_n,yp_n, r, w) { | |
# ## D2 y + 2 r w D y + w^2 y = 0 or # f = f(t, y, yp ) | |
ans = -2.0*r*w*yp_n - w*w*y_n | |
ans | |
} | |
### Runge Kutta 2order, call every step | |
## incomplete! | |
## incomplete! | |
## incomplete! | |
integrate_f <- function(t_n,y_n,yp_n,h=0.0009765, r=1.0, w = 1.0) { | |
one_sixth = 1.0/6.0 | |
half_h = 0.5 *h | |
k1 = h*f(t_n, y_n, yp_n, r, w) | |
k2 = h*f(t_n + half_h, y_n + half_h*y_p, yp_n + 0.5*k1, r, w) | |
k3 = h*f(t_n + half_h, y_n + half_h*(y_p + 0.5*half_h*k1), yp_n + 0.5*k2, r, w) | |
k4 = h* f(t_n + h, y_n + yp_n*h + half_h*k2, yp_n + k3) | |
y_n_new = y_p + h*(yp_n + one_sixth*(k1+k2+k3)) | |
yp_n_new = yp_n + one_sixth*(k1+2*k2+2*k3+k4) | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
##declare equation F(t,x,v)= -(k/m)v-(b/m)x | |
input k,b,m | |
x(0) <- 0 ##declare initial | |
v(0) <- 0 ##declare initial velocity | |
h <- 0.2 ##declare time step | |
ti <- 0 ## declare tinitial | |
tf <- 5 ## declare tfinal | |
x <- x(0) | |
v <- v(0) | |
loop for t from ti to tf with increment of h | |
x <-array(, dim = c(25,2)) | |
dx1 = h*v | |
dv1 = h*F(t,x,v) | |
dx2 = h*(v+dv1/2) | |
dv2 = h*F(t+h/2, x+dx1/2, v+dv1/2) | |
dx3 = h*(v+dv2/2) | |
dv3 = h*F(t+h/2, x+dx2/2, v+dv2/2) | |
dx4 = h*(v+dv3) | |
dv4 = h*F(t+h x+dx3, v+dv3) | |
dx = (dx1 + 2*dx2 + 2*dx3 + dx4)/6 | |
dv = (dv1 + 2*dv2 + 2*dv3 + dv4)/6 | |
x = x+dx | |
v = v+dv | |
end loop |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment