/* 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.mac, Ch.3 Computational Physics with Maxima or R: Boundary Value and Eigenvalue Problems Ted Woollett http://www.csulb.edu/~woollett August, 2015 */ /*** from cp1code.mac: simple search with step halving rtsearch(...) could be improved by counting the number of steps and step-halving, and imposing an upper limit on either or both. Better code is in the code files for Example 3, FW.mac and LJ6-12.mac, and is called bracket_basic and bracket. The output should be combined with find_root to find accurate roots. ****/ rtsearch(func,xx,dxx,xacc) := block([fold, fnew, x:xx, dx:dxx], do ( if abs(dx) <= xacc then return(), fold : func(x), x : x + dx, fnew : func(x), if fold*fnew < 0 then ( x : x - dx, dx : dx/2) ), x)$ /* example of use: (%i1) load(cp3); (%o1) "c:/k3/cp3.mac" (%i2) rtsearch(lambda([x],exp(x)*log(x) - x^2),1.5,0.1,0.01); (%o2) 1.6875 (%i3) rtsearch(lambda([x],exp(x)*log(x) - x^2),1.5,0.1,0.001); (%o3) 1.69375 */ /******* Secant Method search for root, from cp1code.mac and cp1.pdf, Sec. 1.6.9 */ secant(func, vv0, vv1,[oa]) := block([eps0:1e-5,x0,x1 ,x2,ss,jj:0 ,jjmax :3000], if length(oa) > 0 then eps0 : oa[1], x0 : float(vv0), x1 : float(vv1), do (jj : jj + 1, if jj > jjmax then ( print(" exceeded jjmax limit "), return(x1)), if abs(func(x1)) < eps0 then return(x1), ss : func(x1) - func(x0), if equal(ss,0) then ( print(" denominator near",(x0+x1)/2,"is zero "), return(x1) ), x2 : x1 - func(x1)*(x1 - x0)/ss, x0 : x1, x1 : x2))$ /* examples of use (%i4) secant(lambda([x],exp(x)*log(x) - x^2),1.5,1.6,0.001); (%o4) 1.694558 (%i5) secant(sin,2,3,1e-8); (%o5) 3.141593 (%i6) secant(lambda([x],x^2-5),1,2); (%o6) 2.236068 */ /* sho(h,y0,yp0) integrates 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: the list: [ [0, y0], [h, y1], [2*h, y2],...] */ sho(h,y0,yp0) := block([A,y1,N,ym,yz,yp,rL,x, xmin:0,xmax:1,numer],numer:true, A : 2*(1 - 5*%pi^2*h^2/3)/(1 + %pi^2*h^2/3), /* print(" A = ",A), */ y1 : y0 + h*yp0 - 2*%pi^2*h^2*y0 - 2*%pi^2*h^3*yp0/3 + 2*%pi^4*h^4*y0/3, /* print(" y1 = ",y1), */ N : round( (xmax - xmin)/h ), rL : [[xmin+h,y1], [xmin,y0]], x : xmin + 2*h, ym : y0, yz : y1, for j thru N - 1 do ( yp : A*yz - ym, rL : cons ( [x,yp], rL), ym : yz, yz : yp, x : x + h), reverse(rL))$ /* (%i1) load(cp3); (%o1) "c:/k3/cp3.mac" (%i2) soln : sho(0.01, 1, 0)$ (%i3) xL : take(soln,1)$ (%i4) fll(xL); (%o4) [0,1.0,101] (%i5) head(xL); (%o5) [0,0.01,0.02,0.03,0.04,0.05] (%i6) tail(xL); (%o6) [0.95,0.96,0.97,0.98,0.99,1.0] (%i7) xL[1]; (%o7) 0 (%i8) yL : take(soln,2)$ (%i9) fll(yL); (%o9) [1,1.0,101] (%i10) head(yL); (%o10) [1,0.998027,0.992115,0.982287,0.968583,0.951057] (%i11) yL[1]; (%o11) 1 (%i12) plot2d([['discrete,xL,yL],cos(2*%pi*x)],['x,0,1],['xlabel,"x"], ['ylabel,""],['legend,"numerov","analytic"], [gnuplot_preamble,"set key bottom"])$ (%i13) plot2d([['discrete,soln],cos(2*%pi*x)],['x,0,1],['xlabel,"x"], ['ylabel,""],['legend,"numerov","analytic"], [gnuplot_preamble,"set key bottom"])$ */ /* sho2(h,y0,yp0) is designed to integrate the ode of the classical simple harmonic oscillator from x = 0 to x = 1 (like sho(...) above), using the Numerov method discussed in example3.pdf, except that this version uses hashed arrays xL[j] and yL[j] to accumulate the solution. We use the line local(xL,yL) at the top of the program so that we can reuse the names globally (see example below). We use the listarray function to return the list: [ xlist, ylist]. */ sho2(h,y0,yp0) := block([A,y1,xmin,xmax,N,x,ym,yz,yp, numer], numer:true, local(xL,yL), 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, print(" A = ",A," y1 = ", y1), xmin : 0, xmax : 1, N : round ( (xmax - xmin)/h), print(" N = ", N), xL[1] : xmin, xL[2] : xmin + h, yL[1] : y0, yL[2] : y1, x : xmin + 2*h, ym : y0, yz : y1, for j : 3 thru N + 1 do ( yp : A*yz - ym, xL[j] : x, yL[j] : yp, ym : yz, yz : yp, x : x + h ), [ listarray(xL), listarray(yL) ] )$ /* (%i14) remvalue(xL, yL); (%o14) [xL,yL] (%i15) xL[1]; (%o15) xL[1] (%i16) yL[1]; (%o16) yL[1] (%i17) soln : sho2(0.01, 1, 0)$ A = 1.9960535 y1 = 0.998027 N = 100 (%i18) xL[1]; (%o18) xL[1] (%i19) xL : soln[1]$ (%i20) fll(xL); (%o20) [0,1.0,101] (%i21) head(xL); (%o21) [0,0.01,0.02,0.03,0.04,0.05] (%i22) yL[1]; (%o22) yL[1] (%i23) yL : soln[2]$ (%i24) fll(yL); (%o24) [1,1.0,101] (%i25) plot2d([['discrete,xL,yL],cos(2*%pi*x)],['x,0,1],['xlabel,"x"], ['ylabel,""],['legend,"numerov","analytic"], [gnuplot_preamble,"set key bottom"])$ */ /* rk4_step(ode,var,init,domain) if independent variable is t, then domain = [t, t0, dt] to advance the solution one Runge-Kutta step from t0 to (t0+dt). 1D example: to advance the solution of the single first order ode dy/dt = cos(t), assuming that y(0) = 0, from t=0 to t = 0.01 with a single Runge-Kutta step: (%i4) load(cp3); (%o4) "c:/k3/cp3.mac" (%i5) soln : rk4_step(cos(t),y,0,[t,0,0.01]); (%o5) [0.01,0.00999983] (%i6) sin(0.01); (%o6) 0.00999983 2D example: to advance the solution of the pair of first order ode's: dx/dt = vx, d(vx)/dt = -4*%pi^2*x, with the initial conditions: x(0) = 1, vx(0) = 0, and the analytic solution x = cos(2*%pi*t), vx = -2*%pi*sin(2*%pi*t) from t = 0 to t = 0.01 via a single Runge-Kutta step: (%i1) load(cp3); (%o1) "c:/k3/cp3.mac" (%i2) soln : rk4_step([vx,-4*%pi^2*x],[x,vx],[1, 0],[t, 0, 0.01]); (%o2) [0.01,0.998027,-0.394524] (%i3) cos(2*%pi*0.01),numer; (%o3) 0.998027 (%i4) -2*%pi*sin(2*%pi*0.01),numer; (%o4) -0.394524 */ rk4_step(ode, var, init, domain) := block([uvw,k1,k2,k3,k4,t0,t1,dt, numer:true,display2d:false], init : float(init), domain : float(domain), if (not(listp(ode))) then ( ode : [ode], var : [var], init : [init]), local(rkfunc), define(funmake(rkfunc,cons(domain[1],var)),float(ode)), dt : domain[3], t0 : domain[2], t1 : t0 + dt, uvw: init, k1: apply(rkfunc,cons(t0,uvw)), k2: apply(rkfunc,cons((t0+t1)/2, uvw+k1*dt/2)), k3: apply(rkfunc,cons((t0+t1)/2,uvw+k2*dt/2)), k4: apply(rkfunc,cons(t1,uvw+k3*dt)), /* print(" k1 = ",k1," k2 = ",k2), print(" k3 = ",k3," k4 = ",k4), */ uvw: uvw + dt*(k1+2*k2+2*k3+k4)/6, uvw : expand(uvw), cons(t1,uvw))$ /* (%i24) rk4_step(-t*y,y,2,[t,0,0.1]); (%o24) [0.1,1.990025] (%i25) rk4_step([vx,-4*%pi^2*x],[x,vx],[1,0],[t,0,0.1]); (%o25) [0.1,0.8091,-3.688084] */ /* rk4 syntax: similar to rk4_step, except domain = [t, t0, tf, dt] for the same two examples as above, rk4(-t*y,y,2,[t, 0, 1, 0.01]) rk4([vx,-4*%pi^2*x],[x,vx],[1,0],[t, 0, 1, 0.01]); etc. output is the list (for the case of two first order odes, where x1(t) = x(t) and x2(t) = x'(t) = vx(t)), [ [t0,x1(t0),x2(t0)], [t0+h,x1(t0+h),x2(t0+h)],...] See ch. 2 for examples of use. */ rk4(ode, var, init, domain) := block([uvw,rksoln,n,k1,k2,k3,k4,t0,t1,dt, r,numer:true,display2d:false], init : float(init), domain : float(domain), if (not(listp(ode))) then ( ode : [ode], var : [var], init : [init]), /***** added mar 30, 2015 ***/ for i thru length(init) do if not floatnump( init[i] ) then ( print ( " init[",i,"] is not a floating point number"), error("error return from rk4")), /***** ****/ local(rkfunc), define(funmake(rkfunc,cons(domain[1],var)),float(ode)), translate(rkfunc), dt : domain[4], t0 : domain[2], /* n: fix((domain[3] - t0)/dt), */ n: round((domain[3] - t0)/dt), uvw: init, if (not(numberp(last(apply(rkfunc,cons(t0,uvw)))))) then error("Expecting a number when the initial state is replaced in the equations, but instead found:" ,last(apply(rkfunc,cons(t0,uvw)))), rksoln: [cons(t0, init)], for i thru n do ( r: errcatch ( t1: domain[2]+i*dt, k1: apply(rkfunc,cons(t0,uvw)), k2: apply(rkfunc,cons((t0+t1)/2, uvw+k1*dt/2)), k3: apply(rkfunc,cons((t0+t1)/2,uvw+k2*dt/2)), k4: apply(rkfunc,cons(t1,uvw+k3*dt))), if length(r)=0 then return() else uvw: uvw + dt*(k1+2*k2+2*k3+k4)/6, t0: t1, rksoln : cons(cons(t0,uvw), rksoln)), reverse(rksoln))$ /* linear1 ([dy1dx,dy2dx], [y1,y2], [A,B], [x,a,b,h]) 2nd order linear ode BVP method see Ch. 3 of Computational Physics with Maxima or R. 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 x is the independent variable, h is the grid spacing. Two guesses of y'(a) are asked for interactively. Returns a Maxima list with elements of the form [xi, y1i, y2i], which are approximations to [x, y(x), y'(x)] at the grid points. */ linear1(dy1dy2, ivar, bc, xgrid) := block([A,B,del, xL,uL,upL, ub, vL, vpL,vb, outL, sdiff, alpha, y1L, y2L, small : 1e-12, numer],numer:true, if length(dy1dy2) # 2 then ( print(" need two first order ode derivatives"), return(done)), A : bc[1], B : bc[2], /* print(" A = ",A," B = ",B), print(" domain is", xgrid[2], " <=", xgrid[1]," <= ", xgrid[3]), print(" grid interval is ", xgrid[4]), */ /* get first guess for y'(a) */ del : read( " first guess y'(a) = " ), /* print(" del = ",del), */ /* get first IVP solution */ outL : rk4(dy1dy2,ivar,[A,del],xgrid), xL : take(outL,1), uL : take(outL,2), upL : take(outL,3), /* get second guess for y'(a) */ del : read( " second guess y'(a) = " ), /* get second IVP solution */ outL : rk4(dy1dy2,ivar,[A,del],xgrid), vL : take(outL,2), vpL : take(outL,3), ub : last(uL), vb : last(vL), sdiff : ub - vb, /* check that sdiff is not "zero" so we can divide by it */ if abs(sdiff) < small then ( print(" sdiff = ub - vb is too small; choose different values for y'(a) guesses. "), return(done)), /* 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, makelist( [xL[j], y1L[j], y2L[j] ], j, 1, length(xL)))$ /* (%i1) load(cp3); (%o1) "c:/k3/cp3.mac" (%i14) soln : linear1( [y2, -(1+y1)*%pi^2/4], [y1,y2], [0,1], [x,0,1,0.01] )$ first guess y'(a) = 1; second guess y'(a) = 1.5; (%i18) fll(soln); (%o18) [[0.0,0.0,3.141593],[1.0,1.0,-1.570796],101] (%i15) xL : take(soln,1)$ (%i16) yL : take(soln,2)$ (%i17) plot2d([[discrete,xL,yL],cos(%pi*x/2)+2*sin(%pi*x/2) -1], [x,0,1], [legend,"numerical","exact"],[ylabel,"y"], [gnuplot_preamble,"set grid"])$ */ /* shoot1: uses secant search method. see Ch. 3 of Computational Physics with Maxima or R. example of method that works both for linear and nonlinear ODE's; this example is specific to the linear ode: u'' = -(1+u)*pi^2/4, u(0) = 0, u(1) = 1; y1=u, y2=u' . seek solution such that F(del) = y1(del;x=1) -1 = 0 within tolerance tol . shoot1 returns the final Runge-Kutta solution list returned by rk4 */ shoot1(del0,del1,h,tol) := block([del2, F0, F1, F2, outL,k,kmax:10, dy1dy2,ivar,xrange, numer],numer:true, dy1dy2 : [y2, -(1+y1)*%pi^2/4], ivar : [y1,y2], xrange : [x,0,1,h], /* get first trial solution */ outL : rk4(dy1dy2,ivar,[0,del0],xrange), F0 : second(last(outL)) - 1, /* get second trial solution */ outL : rk4(dy1dy2,ivar,[0,del1],xrange), F1 : second(last(outL)) - 1, /* update adjustable parameter del using secant rule */ for k thru kmax do ( del2 : del1 - F1*(del1 - del0)/(F1-F0), outL : rk4(dy1dy2,ivar,[0,del2],xrange), F2 : second(last(outL)) - 1, print(" k = ",k, " del2 = ",del2," F2 = ",F2), if abs(F2) < tol then return(), F0 : F1, /* rotate values */ F1 : F2, del0 : del1, del1 : del2), outL)$ /* (%i3) soln : shoot1(1,1.01,0.01,1e-6)$ k = 1 del2 = 3.141593 F2 = 5.9507954E-14 (%i4) fll(soln); (%o4) [[0.0,0.0,3.141593],[1.0,1.0,-1.570796],101] (%i5) xL : take(soln,1)$ (%i6) fll(xL); (%o6) [0.0,1.0,101] (%i7) y1L : take(soln,2)$ (%i8) fll(y1L); (%o8) [0.0,1.0,101] */ /* shoot2: uses secant search method see Ch. 3 of Computational Physics with Maxima or R. example of method that works both for linear and nonlinear ODE's; 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: y1 = y, y2 = y' . The grid size h = 0.01 is hardwired. 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 Runge-Kutta solution list returned by our homemade rk4. The analytic solution is y(x) = 1/(3+x). */ shoot2(del0,del1,tol) := block([del2, F0, F1, F2, outL, j, jmax:10, h:0.01, yinit:1/2, yfinal:1/3, dy1dy2, ivar, xrange, numer],numer:true, dy1dy2 : [y2, 2*y1^3], ivar : [y1,y2], xrange : [x, -1, 0, h], /* get first trial solution */ outL : rk4(dy1dy2,ivar,[yinit, del0], xrange), F0 : second(last(outL)) - yfinal, /* get second trial solution */ outL : rk4(dy1dy2,ivar,[yinit, del1], xrange), F1 : second(last(outL)) - yfinal, /* update adjustable parameter del using secant rule */ for j thru jmax do ( del2 : del1 - F1*(del1 - del0)/(F1-F0), outL : rk4(dy1dy2,ivar,[yinit, del2], xrange), F2 : second(last(outL)) - yfinal, print(" j = ",j, " del2 = ",del2," F2 = ",F2), if abs(F2) < tol then return(), /* escape from loop */ F0 : F1, /* rotate values */ F1 : F2, del0 : del1, del1 : del2), outL)$ /* (%i3) soln : shoot2(0,-0.5,1e-6)$ j = 1 del2 = -0.26225 F2 = -0.01433 j = 2 del2 = -0.24948 F2 = 6.084312E-4 j = 3 del2 = -0.25 F2 = -1.4480274E-6 j = 4 del2 = -0.25 F2 = -1.4686691E-10 (%i4) fll(soln); (%o4) [[-1.0,0.5,-0.25],[0.0,0.33333,-0.11111],101] (%i5) xL : take(soln,1)$ (%i6) fll(xL); (%o6) [-1.0,0.0,101] (%i7) y1L : take(soln,2)$ (%i8) fll(y1L); (%o8) [0.5,0.33333,101] */ /* simpson's 1/3 rule quadrature for case of discrete data. modified from Ch. 1 code file cp1code.mac */ simp(xv, yv) := block([n, s], n : length(xv) - 1, /* number of steps */ if mod(n,2) # 0 then ( print(" in simp(xv, yv), length of xv and yv should be odd" ), return()), h : xv[2] - xv[1], if n = 2 then s : yv[1] + 4*yv[2] + yv[3] else s : yv[1] + yv[n+1] + 4*apply("+",makelist(yv[i],i,2,n,2)) + 2*apply("+", makelist(yv[i],i,3,n-1,2)), float(s*h/3))$ /* (%i41) xL : makelist(x,x,0,1,0.01)$ (%i42) yL : map('sin,xL)$ (%i43) fll(xL); (%o43) [0,1.0,101] (%i44) fll(yL); (%o44) [0,0.84147,101] (%i45) simp(xL,yL); (%o45) 0.4597 (%i47) integrate(sin(x),x,0,1),numer; (%o47) 0.4597 (%i48) simp( reverse(xL), reverse (yL) ); (%o48) -0.4597 (%i49) integrate(sin(x)^2,x,0,1),numer; (%o49) 0.27268 (%i50) simp(xL,yL^2); (%o50) 0.27268 (%i51) simp(reverse(xL),reverse(yL)^2); (%o51) -0.27268 (%i52) simp(reverse(xL),yL^2); (%o52) -0.27268 */ /* sec. 1.5.8 Trapezoidal Rule: Non-uniform grid this code is from our chapter 1 cp1code.mac file. */ trapz(xv,yv) := block([n:length(xv)], dxx : makelist( xv[i] - xv[i-1], i, 2, n), yy : makelist( yv[i-1] + yv[i], i, 2, n), float(dxx . yy/2))$ /* list utilities */ head(mylist) := block([numL,nleft:6], numL : length(mylist), rest (mylist, -(numL - nleft)))$ tail(mylist) := block([numL,nleft:6], numL : length(mylist), rest (mylist, numL - nleft))$ /* a better name for take would be get_vec */ take(%aL,%nn) := (map(lambda([x],part(x,%nn)), %aL))$ range(aaL) := print(" min = ",lmin(aaL)," max = ", lmax(aaL))$ fll(x) := [first(x),last(x),length(x)]$ first_nonzero(aL) := block([np:0, numer],numer:true, for j thru length(aL) do ( if is(notequal(0.0,aL[j])) then ( np : j, return())), np)$ merge(list1,list2) := (flatten ( cons(list1, list2) ) )$ /* (%i21) aa : [a1,a2,a3]$ (%i22) bb : [b1,b2,b3]$ (%i25) merge(aa,bb); (%o25) [a1,a2,a3,b1,b2,b3] */ ratprint:false$ fpprintprec:8$