## COPYRIGHT NOTICE ## Copyright (C) 2015 Edwin L. (Ted) Woollett ## http://www.csulb.edu/~woollett ## This program is free software; you can redistribute ## it and/or modify it under the terms of the ## GNU General Public License as published by ## the Free Software Foundation; either version 2 ## of the License, or (at your option) any later version. ## This program is distributed in the hope that it ## will be useful, but WITHOUT ANY WARRANTY; ## without even the implied warranty of MERCHANTABILITY ## or FITNESS FOR A PARTICULAR PURPOSE. See the ## GNU General Public License for more details at ## http://www.gnu.org/copyleft/gpl.html ## cp3.R utilities for Ch. 3 of Computational Physics ## with Maxima or R, Boundary Value and Eigenvalue Problems ## myrk4( init, grid, func) ## for a single first order ode (one dependent variable), myrk4 returns a R vector containing ## the values of that dependent variable (and func should return a single derivative value). ## For two or more first order odes (two or more dependent variables), myrk4 returns ## a R list of vectors, one for each dependent variable (and func should return a ## R vector of derivative values). See 1D and 2D examples after the code. ## this code is copied from myode.R from chapter 2. myrk4 = function(init, grid ,func) { num.var = length(init) solnList = list() n.grid = length(grid) for (k in 1:num.var) { solnList[[k]] = vector(length = n.grid) solnList[[k]][1] = init[k] } h = grid[2] - grid[1] # step size yL = vector(length = num.var) # solution at beginning of each step for (j in 1:(n.grid-1)) { for (k in 1:num.var) yL[k] = solnList[[k]][j] k1 = func(grid[j], yL) # vector of derivatives k2 = func(grid[j] + h/2, yL + h*k1/2) # vector of derivatives k3 = func(grid[j] + h/2, yL + h*k2/2) # vector of derivatives k4 = func(grid[j] + h, yL + h*k3) # vector of derivatives for (k in 1:num.var) solnList[[k]][j+1] = solnList[[k]][j] + h*(k1[k] + 2*k2[k] + 2*k3[k] + k4[k])/6 } if (num.var==1) solnList[[1]] else solnList} ## two dependent variables example of use of myrk4. ## > source("cp3.R") ## > tL = seq(0,1,0.1) ## > fll(tL) ## 0 1 11 ## note that the func sho returns a R vector of derivative values. ## > sho = function(t,y) { c( y[2], -4*pi^2*y[1] ) } ## > out = myrk4(c(1,0), tL, sho) ## > xL = out[[1]] ## > xL ## [1] 1.000000 0.809102 0.310104 -0.306633 -0.806047 -0.997964 -0.809517 -0.312810 0.302669 ## [10] 0.802336 0.995920 ## > vxL = out[[2]] ## > vxL ## [1] 0.0000000 -3.6880842 -5.9680715 -5.9724674 -3.7014458 -0.0220779 3.6627122 5.9490744 ## [9] 5.9670776 3.7117057 0.0440659 ## > tL ## [1] 0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 ## one dependent variable example using myrk4: ## dy/dt = -y, y(0)=1 solution vector using function deriv and same tL grid of times. ## > deriv = function(t, y) {-y} ## > out = myrk4(1,tL,deriv) ## > is.vector(out) ## [1] TRUE ## > out ## [1] 1.000000 0.904837 0.818731 0.740818 0.670320 0.606531 0.548812 0.496586 0.449329 0.406570 ## [11] 0.367880 ## ## ====================================================== ## linear1 (bc, grid, func): see Ch. 3 of Computational Physics ## with Maxima or R discussion. ## linear 2nd order ode BVP method over finite interval ## for case: ## y''(x) = p(x) y'(x) + q(x) y(x) + r(x), a <= x <= b ## y1(x) = y(x), y2(x) = y'(x) ## y(a) = A, y(b) = B ## Two guesses of y'(a) are asked for interactively. ## The vector 'grid' can be generated as xL = seq(a,b,h), for example. ## linear1 calls homemade myrk4, which expects func to return ## a R vector c(y[2], y'' as a function of x,y[1], and y[2]). ## For the above problem, bc = c(A,B). ## linear1 returns list(y1L, y2L), in which y1L is a R vector of y(x) grid ## values, and y2L is a R vector of y'(x) grid values. linear1 = function(bc, grid, func) { A = bc[1] B = bc[2] ## get first guess for y'(a) del = as.numeric ( readline ( " first guess y'(a) = " )) ## get first IVP solution ; myrk4 returns list(y1L, y2L) outL = myrk4 ( init = c(A, del), grid = grid, func = func ) uL = outL[[1]] upL = outL[[2]] ## get second guess for y'(a) del = as.numeric ( readline ( " second guess y'(a) = " )) ## get second IVP solution outL = myrk4 ( init = c(A, del), grid = grid, func = func ) vL = outL[[1]] vpL = outL[[2]] ub = last(uL) vb = last(vL) sdiff = ub - vb ## check that sdiff is not "zero" so we can divide by it if (abs(sdiff) < 1e-12) { cat(" sdiff = ub - vb is too small; choose different values for y'(a) guesses. \n ") return() } ## form numerical solution y1L and its first derivative y2L alpha = (B - vb)/sdiff # scalar number needed to construct solution y1L = alpha*uL + (1 - alpha)*vL y2L = alpha*upL + (1 - alpha)*vpL list ( y1L, y2L ) } ## use of linear1 ## solves y'' = -pi^2 /4 * (1 + y), for y(0) = 0, y(1) = 1 ## > source("cp3.R") ## > ex1 = function(x,y) { c( y[2], -(1 + y[1])*pi^2/4) } ## > xL = seq(0, 1, 0.01) ## > BC = c(0, 1) ## > soln = linear1(bc=BC, grid = xL, func = ex1) ## first guess y'(a) = 1 ## second guess y'(a) = 1.5 ## > yL = soln[[1]] ## > fll(yL) ## 0 1 101 ## > ypL = soln[[2]] ## > fll(ypL) ## 3.141593 -1.570796 101 ## > plot(xL,yL, type = "l", lwd = 2, col = "blue", ylab = "y(x)", xlab = "x") ## > mygrid() ## > curve ( cos(pi*x/2)+2*sin(pi*x/2) -1,0,1,add = TRUE,col = "red") ## > legend("bottomright", col = c("blue","red"), legend = c( "numerical", "exact"), ## + lwd = 2) ## shoot2: uses secant search method: see Ch. 3 of Computational Physics ## with Maxima or R discussion. ## example of method that works both for linear and nonlinear ## ODE BVPs; this example is specific to the nonlinear ode: ## y'' = 2*y^3, -1 <= x <= 0, ## with boundary conditions: y(-1) = 1/2, y(0) = 1/3. ## Use variables: y[1] = y, y[2] = y' . ## Uses secant search method to ## seek a value of del = y'(-1) such that F(del) = y1(del, x=0) -1/3 = 0 within tolerance tol . ## shoot2 returns the final solution list(yL, ypL) returned by our homemade myrk4. ## The analytic solution is y(x) = 1/(3+x). ## Needs predefined derivatives function derivs, and predefined grid values xL. shoot2 = function(del0,del1,tol) { jmax = 10 yinit = 1/2 yfinal = 1/3 ## get first trial solution outL = myrk4(init = c(yinit,del0), grid = xL, func = derivs) F0 = last(outL[[1]]) - yfinal ## get second trial solution outL = myrk4(init = c(yinit,del1), grid = xL, func = derivs) F1 = last(outL[[1]]) - yfinal ## update adjustable parameter del using secant rule for ( j in seq(1,jmax,1) ) { del2 = del1 - F1*(del1 - del0)/(F1-F0) outL = myrk4(init = c(yinit,del2), grid = xL, func = derivs) F2 = last(outL[[1]]) - yfinal cat(" j = ",j, " del2 = ",del2," F2 = ",F2, "\n") if ( abs(F2) < tol ) break ## escape from loop F0 = F1 ## rotate values F1 = F2 del0 = del1 del1 = del2} outL} ## > source("cp3.R") ## > derivs = function(x,y) { c(y[2], 2*y[1]^3) } ## > xL = seq(-1, 0, 0.01) ## > soln = shoot2(0,-0.5,1e-6) ## j = 1 del2 = -0.2622537 F2 = -0.01432995 ## j = 2 del2 = -0.249481 F2 = 0.0006084312 ## j = 3 del2 = -0.2500012 F2 = -1.448027e-06 ## j = 4 del2 = -0.25 F2 = -1.468672e-10 ## > is.list(soln) ## [1] TRUE ## > y1L = soln[[1]] ## > fll(y1L) ## 0.5 0.3333333 101 ## > plot(xL, y1L,type="l", col = "blue", lwd = 3,xlab = "x", ylab = "y") ## > mygrid() ## > curve(1/(x + 3),-1,0,col = "red", lwd = 3, add = TRUE) ## > legend("topright", col = c("blue", "red"), lwd = 2, ## + legend = c("numerical", "analytic"), cex = 1.5) ## sho(h,y0,yp0): see Ch. 3 of Computational Physics ## with Maxima or R discussion. ## integrates the classical simple harmonic oscillator with unit period ## d^2y/dx^2 = - 4 pi^2 y(x) over [x,0,1] ## using the Numerov method discussed in example3.pdf. ## Input: h = step size, y0 = y(0), yp0 = y'(0) ## Output: list( xL, yL) sho = function(h,y0,yp0) { A = 2*(1 - 5*pi^2*h^2/3)/(1 + pi^2*h^2/3) y1 = y0 + h*yp0 - 2*pi^2*h^2*y0 - 2*pi^2*h^3*yp0/3 + 2*pi^4*h^4*y0/3 cat (" A = ", A, " y1 = ", y1, "\n") xmin = 0 xmax = 1 N = round( (xmax - xmin)/h ) cat (" N = ", N, "\n") yL = vector(length = N + 1) xL = vector(length = N + 1) xL[1] = xmin xL[2] = xmin + h yL[1] = y0 yL[2] = y1 x = xmin + 2*h ym = y0 yz = y1 for (j in 3:(N + 1)) { yp = A*yz - ym xL[j] = x yL[j] = yp ym = yz yz = yp x = x + h } list( xL, yL) } ## > source("cp3.R") ## > soln = sho(0.01,1,0) ## A = 1.99605 y1 = 0.998027 ## N = 100 ## > xL[1] ## Error: object 'xL' not found ## > xL = soln[[1]] ## > fll(xL) ## 0 1 101 ## > yL = soln[[2]] ## > fll(yL) ## 1 1 101 ## rk4_step(init,t0,dt,func) returns y(t0+dt) if one dependent variable ## if one dependent variable y(t) , init = y(t0), t0 = initial value of ## independent variable, dt = Runge-Kutta step size, func = derivative function ## or returns the vector c(x(t+dt), vx(t+dt)) if simple harmonic oscillator (two dependent ## variables) (see example below), etc. rk4_step = function(init, t, dt, func) { num.var = length(init) h = dt # step size yL = init k1 = func(t, yL) # vector of derivatives k2 = func(t + h/2, yL + h*k1/2) # vector of derivatives k3 = func(t + h/2, yL + h*k2/2) # vector of derivatives k4 = func(t + h, yL + h*k3) # vector of derivatives for (k in 1:num.var) yL[k] = yL[k] + h*(k1[k] + 2*k2[k] + 2*k3[k] + k4[k])/6 if (num.var==1) yL[1] else yL} ## example 1: one dependent variable, integrate dy/dt = -y, one step dt = 0.01 ## from t=0, y(0)=1, analytic soln is y(t) = exp(-t). ## > deriv = function(t, y) {-y} ## > rk4_step(1,0,0.01,deriv) ## [1] 0.99005 ## > exp(-0.01) ## [1] 0.99005 ## example2: two dependent variables, simple harmonic oscillator with unit period ## dx/dt=vx = y[2], d(vx)/dt = -4*pi^2*x, x = y[1], with t0 = 0, dt = 0.01, x(0)=1, vx(0) = 0. ## init = c(x(t0), vx(t0)) = c( y10, y20 ). ## > sho = function(t,y) { c( y[2], -4*pi^2*y[1] ) } ## > rk4_step(c(1,0),0,0.01,sho) ## [1] 0.998027 -0.394524 ## > cos(2*pi*0.01) ## [1] 0.998027 ## > -2*pi*sin(2*pi*0.01) ## [1] -0.394524 ## Trapezoidal Rule: allows non-uniform grid, from Ch. 1 file cp1code.R trapz = function(xv,yv) { idx = 2:length(xv) as.numeric( (xv[idx] - xv[idx - 1]) %*% (yv[idx - 1] + yv[idx])/2)} ## modified from cp1code.R: Simpson's 1/3 Rule ## simp(xv,yv) for discrete data vectors simp = function(xv,yv) { if ( length(xv) != length(yv) ) { cat(" in simp(xv, yv) vectors xv and yv must have the same length \n") return() } N = length(xv) - 1 if (N %% 2 != 0) { cat(" in simp(xv,yv) both xv and yv should have length = odd \n") return()} h = xv[2] - xv[1] if (N == 2) s = yv[1]+4*yv[2]+yv[3] else { s = yv[1] + yv[N+1] + 4*sum(yv[seq(2,N,by=2)]) + 2*sum(yv[seq(3,N-1,by=2)])} s*h/3} ## simp2(func,x1,x2,N) numerical integration integrate(func(x),x1,x2) ## using N panels; N should be even integer simp2 = function(func, a, b, N) { if (N %% 2 != 0) { cat(" in simp2(func,a,b,N), N should be an even integer \n") return() } h = (b - a)/N xv = seq(a, b, by=h) yv = func(xv) if (N == 2) s = yv[1]+4*yv[2]+yv[3] else { s = yv[1] + yv[N+1] + 4*sum(yv[seq(2,N,by=2)]) + 2*sum(yv[seq(3,N-1,by=2)])} s*h/3} ## > simp2(sin,0,pi,100) ## [1] 2 ## fll is a home-made vector utility: prints first, last, and length of a vector fll = function(xL) { xlen = length(xL) cat(" ",xL[1]," ",xL[xlen]," ",xlen,"\n") } ## > fll(xL) ## 1 0.9959199 11 mygrid = function() grid(lty = "solid", col = "darkgray") last = function(xL) { xL[ length(xL)] } ## get_last(list.of.vecs) returns a vector whose elements ## are the last components of the vectors get_last = function(rkvecs) sapply(rkvecs, last) ## > myL = list(c(1,2,3,4),c(5,6,7,8),c(9,10,11,12)) ## > get_last(myL) ## [1] 4 8 12 ## rest(xL) returns a vector without the first element rest = function(x) { x[2:length(x)] } ## > aL ## [1] 1 2 3 ## > rest(aL) ## [1] 2 3 options(digits=6) debug1 = FALSE