c c---------------------------------------------------------------- c c problem 1 - see companion paper c implicit real*8 (a-h,o-z) double precision fspace(2000), zeta(4), tol(2), z(4), u(4), err(4) integer ispace(200), m(1), ipar(11), ltol(2) external fsub, dfsub, gsub, dgsub c write (6,99) 99 format(1h1, 35h example of a simple problem setup. . / 46h uniformly loaded beam of variable stiffness, . / 32h simply supported at both ends. /) c c one differential equation of order 4. m(1) = 4 c give location of boundary conditions zeta(1) = 1.d0 zeta(2) = 1.d0 zeta(3) = 2.d0 zeta(4) = 2.d0 c set up parameter array. c use default values for all parameters except for initial c mesh size, no. of tolerances and sizes of work arrays do 10 i=1,11 10 ipar(i) = 0 ipar(3) = 1 ipar(4) = 2 ipar(5) = 2000 ipar(6) = 200 c two error tolerances (on u and its second derivative) ltol(1) = 1 ltol(2) = 3 tol(1) = 1.d-7 tol(2) = 1.d-7 c call colsys (1, m, 1.d0, 2.d0, zeta, ipar, ltol, tol, . dummy, ispace, fspace, iflag, fsub, . dfsub, gsub, dgsub, dummy) c if (iflag .ne. 1) stop c calculate the error at 101 points using the known c exact solution x = 1.d0 do 20 i=1,4 20 err(i) = 0.d0 do 40 j=1,101 call appsln (x, z, fspace, ispace) call exact (x, u) do 30 i=1,4 30 err(i) = dmax1(err(i), dabs(u(i)-z(i))) 40 x = x + .01 write(6,100) (err(i),i=1,4) 100 format(/27h error tolerances satisfied//22h the exact errors are, . / 7x,4d12.4) stop end c................................................................ subroutine fsub (x, z, f) double precision z(4), f(1), x f(1) = (1.d0 - 6.d0*x**2*z(4) - 6.d0*x*z(3)) / x**3 return end c................................................................ subroutine dfsub (x, z, df) double precision z(4), df(1,4), x df(1,1) = 0.d0 df(1,2) = 0.d0 df(1,3) = -6.d0/x**2 df(1,4) = -6.d0/x return end c................................................................ subroutine gsub (i, z, g) double precision z(4), g go to (1, 2, 1, 2), i 1 g = z(1) - 0.d0 return 2 g = z(3) - 0.d0 return end c................................................................ subroutine dgsub (i, z, dg) double precision z(4), dg(4) do 10 j=1,4 10 dg(j) = 0.d0 go to (1, 2, 1, 2), i 1 dg(1) = 1.d0 return 2 dg(3) = 1.d0 return end c................................................................ subroutine exact(x, u) implicit real*8 (a-h,o-z) double precision u(4) c exact solution u(1) = .25d0* (10.d0*dlog(2.d0)-3.d0) * (1.d0-x) + . .5d0* (1.d0/x+ (3.d0+x)*dlog(x) - x) u(2) = -.25d0* (10.d0*dlog(2.d0) - 3.d0) + .5d0 * . (-1.d0/x/x + dlog(x) + (3.d0+x)/x - 1.d0) u(3) = .5d0 * (2.d0/x**3 + 1.d0/x -3.d0/x/x) u(4) = .5d0 * (-6.d0/x**4 - 1.d0/x/x + 6.d0/x**3) return end c c---------------------------------------------------------------- c c problem 2 - see companion paper c implicit real*8 (a-h,o-z) double precision zeta(4), fspace(40000), tol(4), z(4) integer m(2), ipar(11), ispace(2500), ltol(4) common eps, dmu, eps4mu, gamma, xt external solutn,fsub,dfsub,gsub,dgsub c define constants, print a heading. gamma = 1.1d0 eps = .01d0 dmu = eps eps4mu = eps**4/dmu xt = dsqrt(2.d0*(gamma-1.d0)/gamma) write (6,100) gamma, xt, eps, dmu, eps4mu 100 format(1h1,27hdimpling of spherical caps. ./8h gamma =,f7.2/ 6h xt =,d12.5/ 6h eps =,d12.5/ 6h mu =,d12.5/ . 12h eps**4/mu =,d12.5) c define no. of differential equations. ncomp = 2 c orders m(1) = 2 m(2) = 2 c interval ends aleft = 0.d0 aright = 1.d0 c locations of side conditions zeta(1) = 0.d0 zeta(2) = 0.d0 zeta(3) = 1.d0 zeta(4) = 1.d0 c ipar values c a nonlinear problem ipar(1) = 1 c 4 collocation points per subinterval ipar(2) = 4 c initial uniform mesh of 10 subintervals ipar(3) = 10 ipar(8) = 0 c dimension of real work array fspace is 40000 ipar(5) = 40000 c dimension of integer work array ispace is 2500 ipar(6) = 2500 c (these dimensions of fspace and ispace c enable colsys to use meshes of up to 192 intervals.) c print full output. ipar(7) = -1 c initial approximation for nonlinear iteration is provided c in solutn ipar(9) = 1 c a regular problem ipar(10) = 0 c no fixed points in the mesh ipar(11) = 0 c tolerances on all components ipar(4) = 4 do 10 i=1,4 ltol(i) = i 10 tol(i) = 1.d-5 c call colsys call colsys (ncomp, m, aleft, aright, zeta, ipar, ltol, . tol, fixpnt, ispace, fspace, iflag, . fsub,dfsub,gsub,dgsub,solutn) c print values of the obtained approximate solution at points c x = 0,.05, ..., 1. x = 0.d0 write (6,201) 201 format (1h1,44h x phi dphi , . 23h psi dpsi /) 202 format (6x, f5.2, 4x, 6d15.5) np1 = 21 do 555 iii=1,np1 call appsln (x,z,fspace,ispace) write (6,202) x, z x = x + .05d0 555 continue stop end c................................................................ subroutine solutn (x, z, dmval) implicit real*8 (a-h,o-z) common eps, dmu, eps4mu, gamma, xt dimension z(4) , dmval(2) cons = gamma * x * (1.d0-.5d0*x*x) dcons = gamma * (1.d0 - 1.5d0*x*x) d2cons = -3.d0 * gamma * x if (x .gt. xt) go to 10 z(1) = 2.d0 * x z(2) = 2.d0 z(3) = -2.d0*x + cons z(4) = -2.d0 + dcons dmval(2) = d2cons go to 20 10 z(1) = 0.d0 z(2) = 0.d0 z(3) = -cons z(4) = -dcons dmval(2) = -d2cons 20 dmval(1) = 0.d0 return end c................................................................ subroutine fsub (x, z, f) implicit real*8 (a-h,o-z) dimension z(4), f(2) common eps, dmu, eps4mu, gamma, xt f(1) = z(1)/x/x - z(2)/x + (z(1) - z(3)*(1.d0-z(1)/x) - . gamma*x*(1.d0-x*x/2.)) / eps4mu f(2) = z(3)/x/x - z(4)/x + z(1)*(1.d0-z(1)/2.d0/x) / dmu return end c................................................................ subroutine dfsub (x, z, df) implicit real*8 (a-h,o-z) dimension z(4), df(2,4) common eps, dmu, eps4mu, gamma, xt df(1,1) = 1.d0/x/x +(1.d0 + z(3)/x) / eps4mu df(1,2) = -1.d0/x df(1,3) = -(1.d0-z(1)/x) / eps4mu df(1,4) = 0.d0 df(2,1) = (1.d0 - z(1)/x) / dmu df(2,2) = 0.d0 df(2,3) = 1.d0/x/x df(2,4) = -1.d0/x return end c................................................................ subroutine gsub (i, z, g) implicit real*8 (a-h,o-z) dimension z(4) go to (1, 2, 1, 3), i 1 g = z(1) return 2 g = z(3) return 3 g = z(4) - .3d0*z(3) + .7d0 return end c................................................................ subroutine dgsub (i, z, dg) implicit real*8 (a-h,o-z) dimension z(4), dg(4) do 10 j=1,4 10 dg(j) = 0. d0 go to (1, 2, 1, 3), i 1 dg(1) = 1.d0 return 2 dg(3) = 1.d0 return 3 dg(4) = 1.d0 dg(3) = -.3d0 return end c c---------------------------------------------------------------- c c problem 3 - see companion paper c implicit real*8 (a-h,o-z) double precision zeta(5), fspace(40000), tol(2), sval(3), elval(3) integer m(2), ispace(2500), ltol(2), ipar(11) double precision z(5) common en, s, el, cons external fsub, dfsub, gsub, dgsub, solutn data sval/.2d0, .1d0, .05d0/, elval/60.d0, 120.d0, 200.d0/ c en = .2d0 cons = .5d0 * (3.d0-en) ncomp = 2 m(1) = 2 m(2) = 3 aleft = 0.d0 aright = 1.d0 c zeta(1) = 0.d0 zeta(2) = 0.d0 zeta(3) = 0.d0 zeta(4) = 1.d0 zeta(5) = 1.d0 c ipar(1) = 1 ipar(2) = 4 ipar(3) = 10 ipar(4) = 2 ipar(5) = 40000 ipar(6) = 2500 ipar(7) = -1 ipar(8) = 0 ipar(9) = 1 ipar(10) = 0 ipar(11) = 0 c ltol(1) = 1 ltol(2) = 3 tol(1) = 1.d-5 tol(2) = 1.d-5 c c solve a chain of 3 problems do 777 ijk = 1,3 s = sval(ijk) el = elval(ijk) if (ijk .eq. 1) go to 701 c set continuation parameters ipar(9) = 3 ipar(3) = ispace(1) 701 continue write (6,100) en, s, el 100 format(1h1,38h rotating flow over a stationary disk. . /19h parameters - n =,f5.2, 6h s =,f5.2, . 6h l =,f6.1/) c call colsys (ncomp, m, aleft, aright, zeta, ipar, ltol, . tol, fixpnt, ispace, fspace, iflag, . fsub, dfsub, gsub, dgsub, solutn) c if (iflag .ne. 1) stop c print values of the obtained approximate solution at points c x = 0,1,2, ..., l. is6 = ispace(6) is5 = ispace(1) + 2 x = 0.d0 write (6,201) 201 format (1h1,44h x g dg , . 38h h dh d2h/) 202 format (6d15.5) np1 = el + 1.5d0 do 555 iii=1,np1 call approx (ii,x,z,fspace(is6),fspace(1),ispace(1), . fspace(is5),ispace(2),ncomp,m,ispace(4),1,dm,0) xl = x * el z(2) = z(2) / el z(4) = z(4) / el z(5) = z(5) / el / el write (6,202) xl, z x = x + 1.d0 / el 555 continue 777 continue stop end c................................................................ subroutine solutn (x, z, dmval) implicit real*8 (a-h,o-z) common en, s, el, cons double precision z(5), dmval(2) ex = dexp(-el*x) z(1) = 1.d0 - ex z(2) = el * ex z(3) = -el**2 * x**2 * ex z(4) = (el**3 *x**2 - 2.d0 * el**2 * x) * ex z(5) = (-el**4 * x**2 + 4.d0 * el**3 * x - 2.d0 * el**2)*ex dmval(1) = -el * z(2) dmval(2) = (el**5*x*x - 6.d0*el**4*x + 6.d0*el**3) * ex return end c................................................................ subroutine fsub (x, z, f) implicit real*8 (a-h,o-z) double precision z(1), f(1) common en, s, el, cons f(1) = -el * (cons * z(3) * z(2) + (en-1.d0) * z(4) * z(1)) . + el**2 * s * (z(1)-1.d0) f(2) = -el * (cons * z(3) * z(5) + en * z(4)**2) + . el**2 * s * z(4) + el**3 * (1.d0-z(1)**2) return end c................................................................ subroutine dfsub (x, z, df) implicit real*8 (a-h,o-z) double precision z(1), df(2,1) common en, s, el, cons df(1,1) = -el * (en-1.d0) * z(4) + el**2 * s df(1,2) = -el * cons * z(3) df(1,3) = -el * cons * z(2) df(1,4) = -el * (en-1.d0) * z(1) df(1,5) = 0.d0 df(2,1) = -el**3 * 2.d0 * z(1) df(2,2) = 0.d0 df(2,3) = -el * cons * z(5) df(2,4) = -el * en * 2.d0 * z(4) + el**2 * s df(2,5) = -el * cons * z(3) return end c................................................................ subroutine gsub (i, z, g) double precision z(1), g go to (1, 2, 3, 4, 3), i 1 g = z(1) return 2 g = z(3) return 3 g = z(4) return 4 g = z(1) - 1.d0 return end c................................................................ subroutine dgsub (i, z, dg) double precision z(1), dg(1) do 10 j=1,5 10 dg(j) = 0.d0 go to (1, 2, 3, 1, 3), i 1 dg(1) = 1.d0 return 2 dg(3) = 1.d0 return 3 dg(4) = 1.d0 return end