c c---------------------------------------------------------------- c c problem 1 - see companion paper c real 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,99999) c c one differential equation of order 4. m(1) = 4 c give location of boundary conditions zeta(1) = 1. zeta(2) = 1. zeta(3) = 2. zeta(4) = 2. 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 ipar(i) = 0 10 continue 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.e-7 tol(2) = 1.e-7 c call colsys(1, m, 1., 2., 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. do 20 i=1,4 err(i) = 0. 20 continue do 40 j=1,101 call appsln(x, z, fspace, ispace) call exact(x, u) do 30 i=1,4 err(i) = amax1(err(i),abs(u(i)-z(i))) 30 continue x = x + .01 40 continue write (6,99998) (err(i),i=1,4) stop 99999 format (1h1, 35h example of a simple problem setup./10h uniforml, * 36hy loaded beam of variable stiffness,/21h simply supported at, * 11h both ends./) 99998 format (/27h error tolerances satisfied//22h the exact errors are, * /7x, 4e12.4) end c................................................................ subroutine fsub(x, z, f) real z(4), f(1) f(1) = (1.-6.*x**2*z(4)-6.*x*z(3))/x**3 return end c................................................................ subroutine dfsub(x, z, df) real z(4), df(1,4) df(1,1) = 0. df(1,2) = 0. df(1,3) = -6./x**2 df(1,4) = -6./x return end c................................................................ subroutine gsub(i, z, g) real z(4) go to (10, 20, 10, 20), i 10 g = z(1) - 0. return 20 g = z(3) - 0. return end c................................................................ subroutine dgsub(i, z, dg) real z(4), dg(4) do 10 j=1,4 dg(j) = 0. 10 continue go to (20, 30, 20, 30), i 20 dg(1) = 1. return 30 dg(3) = 1. return end c................................................................ subroutine exact(x, u) real u(4) c exact solution u(1) = .25*(10.*alog(2.)-3.)*(1.-x) + .5*(1./x+(3.+x)*alog(x)-x) u(2) = -.25*(10.*alog(2.)-3.) + .5*(-1./x/x+alog(x)+(3.+x)/x-1.) u(3) = .5*(2./x**3+1./x-3./x/x) u(4) = .5*(-6./x**4-1./x/x+6./x**3) return end c c---------------------------------------------------------------- c c problem 2 - see companion paper c real 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.1 eps = .01 dmu = eps eps4mu = eps**4/dmu xt = sqrt(2.*(gamma-1.)/gamma) write (6,99999) gamma, xt, eps, dmu, eps4mu c define no. of differential equations. ncomp = 2 c orders m(1) = 2 m(2) = 2 c interval ends aleft = 0. aright = 1. c locations of side conditions zeta(1) = 0. zeta(2) = 0. zeta(3) = 1. zeta(4) = 1. 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 tol(i) = 1.e-5 10 continue 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. write (6,99998) np1 = 21 do 20 iii=1,np1 call appsln(x, z, fspace, ispace) write (6,99997) x, z x = x + .05 20 continue stop 99999 format (1h1, 27hdimpling of spherical caps./8h gamma =, * f7.2/6h xt =, e12.5/6h eps =, e12.5/6h mu =, e12.5/9h eps**4/m, * 3hu =, e12.5) 99998 format (1h1, 44h x phi dphi , * 23h psi dpsi/) 99997 format (6x, f5.2, 4x, 6e15.5) end c................................................................ subroutine solutn(x, z, dmval) common eps, dmu, eps4mu, gamma, xt dimension z(4), dmval(2) cons = gamma*x*(1.-.5*x*x) dcons = gamma*(1.-1.5*x*x) d2cons = -3.*gamma*x if (x.gt.xt) go to 10 z(1) = 2.*x z(2) = 2. z(3) = -2.*x + cons z(4) = -2. + dcons dmval(2) = d2cons go to 20 10 z(1) = 0. z(2) = 0. z(3) = -cons z(4) = -dcons dmval(2) = -d2cons 20 dmval(1) = 0. return end c................................................................ subroutine fsub(x, z, f) 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.-z(1)/x)-gamma*x*(1.-x*x/ * 2.))/eps4mu f(2) = z(3)/x/x - z(4)/x + z(1)*(1.-z(1)/2./x)/dmu return end c................................................................ subroutine dfsub(x, z, df) dimension z(4), df(2,4) common eps, dmu, eps4mu, gamma, xt df(1,1) = 1./x/x + (1.+z(3)/x)/eps4mu df(1,2) = -1./x df(1,3) = -(1.-z(1)/x)/eps4mu df(1,4) = 0. df(2,1) = (1.-z(1)/x)/dmu df(2,2) = 0. df(2,3) = 1./x/x df(2,4) = -1./x return end c................................................................ subroutine gsub(i, z, g) dimension z(4) go to (10, 20, 10, 30), i 10 g = z(1) return 20 g = z(3) return 30 g = z(4) - .3*z(3) + .7 return end c................................................................ subroutine dgsub(i, z, dg) dimension z(4), dg(4) do 10 j=1,4 dg(j) = 0. 10 continue go to (20, 30, 20, 40), i 20 dg(1) = 1. return 30 dg(3) = 1. return 40 dg(4) = 1. dg(3) = -.3 return end c c---------------------------------------------------------------- c c problem 3 - see companion paper c real zeta(5), fspace(40000), tol(2), sval(3), elval(3) integer m(2), ispace(2500), ltol(2), ipar(11) real z(5) common en, s, el, cons external fsub, dfsub, gsub, dgsub, solutn data sval /.2,.1,.05/, elval /60.,120.,200./ c en = .2 cons = .5*(3.-en) ncomp = 2 m(1) = 2 m(2) = 3 aleft = 0. aright = 1. c zeta(1) = 0. zeta(2) = 0. zeta(3) = 0. zeta(4) = 1. zeta(5) = 1. c ipar(1) = 1 ipar(2) = 4 ipar(3) = 10 ipar(4) = 2 ipar(5) = 40000 ipar(6) = 2500 ipar(7) = 0 ipar(8) = 0 ipar(9) = 1 ipar(10) = 0 ipar(11) = 0 c ltol(1) = 1 ltol(2) = 3 tol(1) = 1.e-5 tol(2) = 1.e-5 c c solve a chain of 3 problems do 30 ijk=1,3 s = sval(ijk) el = elval(ijk) if (ijk.eq.1) go to 10 c set continuation parameters ipar(9) = 3 ipar(3) = ispace(1) 10 continue write (6,99999) en, s, el 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. write (6,99998) np1 = el + 1.5 do 20 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,99997) xl, z x = x + 1./el 20 continue 30 continue stop 99999 format (1h1, 38h rotating flow over a stationary disk./8h parame, * 11hters - n =, f5.2, 6h s =, f5.2, 6h l =, f6.1/) 99998 format (1h1, 44h x g dg , * 38h h dh d2h/) 99997 format (6e15.5) end c................................................................ subroutine solutn(x, z, dmval) common en, s, el, cons real z(5), dmval(2) ex = exp(-el*x) z(1) = 1. - ex z(2) = el*ex z(3) = -el**2*x**2*ex z(4) = (el**3*x**2-2.*el**2*x)*ex z(5) = (-el**4*x**2+4.*el**3*x-2.*el**2)*ex dmval(1) = -el*z(2) dmval(2) = (el**5*x*x-6.*el**4*x+6.*el**3)*ex return end c................................................................ subroutine fsub(x, z, f) real z(1), f(1) common en, s, el, cons f(1) = -el*(cons*z(3)*z(2)+(en-1.)*z(4)*z(1)) + el**2*s*(z(1)-1.) f(2) = -el*(cons*z(3)*z(5)+en*z(4)**2) + el**2*s*z(4) + * el**3*(1.-z(1)**2) return end c................................................................ subroutine dfsub(x, z, df) real z(1), df(2,1) common en, s, el, cons df(1,1) = -el*(en-1.)*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.)*z(1) df(1,5) = 0. df(2,1) = -el**3*2.*z(1) df(2,2) = 0. df(2,3) = -el*cons*z(5) df(2,4) = -el*en*2.*z(4) + el**2*s df(2,5) = -el*cons*z(3) return end c................................................................ subroutine gsub(i, z, g) real z(1), g go to (10, 20, 30, 40, 30), i 10 g = z(1) return 20 g = z(3) return 30 g = z(4) return 40 g = z(1) - 1. return end c................................................................ subroutine dgsub(i, z, dg) real z(1), dg(1) do 10 j=1,5 dg(j) = 0. 10 continue go to (20, 30, 40, 20, 40), i 20 dg(1) = 1. return 30 dg(3) = 1. return 40 dg(4) = 1. return end