sinc-galerkin example

BlockMatrix Sinc-Galerkin Solution ofthe Wind-Driv enCurren tProblem SanoeKoonprasert andKenneth L.Bowers Department ofMathematic alScienc es,Montana StateUniversity Bozeman, MT59717, USA Abstract Ablockmatrix formulation ispresen tedfortheSinc-Galerkin technique applied to thewind-driv encurren tproblem fromoceanograph y.Theblockmatrix formisused todetermine anapproximate solution forthecoupled system ofdieren tialequations describing thewind-driv encurren tproblem. Theapproximate solution determined bythisnewSinc-Galerkin scheme iscompared totheapproximate solution for theoriginal complex-v alued Sinc-Galerkin approac htooceancurren tmodels.The approximate solution isalsocompared toexact solutions toillustrate theexponential convergence rateofthemetho d. Keywords:Blockmatrix, Sinc-Galerkin metho d,wind-driv encurren ts 1Introduction Inrecentyears, thecomplex-v alued Sinc-Galerkin metho dforwind-driv en curren tmodelshasbeenapplied tothemodelformulated asacomplex-v alued ordinary dieren tialequation. However,aSinc-Galerkin technique canalso beapplied toaformulation ofthemodelthatconsists oftworeal-valued coupled ordinary dieren tialequations. Theblockmatrix dependsonthesinc basisfunctions, thedieren tialequations, andtheboundary conditions. The complex-v alued Sinc-Galerkin metho dresults inasmaller discrete system but requires complex arithmetic atallsteps. Theblockformulation eliminates the needforanycomplex-v alued calculations attheexpenseofaslightlylarger matrix system. Inthispaper,weperform anumerical study onamodelofwind-driv encurren ts intheocean. Thismodelisfound in[1].There theSinc-Galerkin technique wasapplied toamodeldescribing wind-driv ensubsurface curren tsincoastal regions andsemi-enclosed seaswhere thevertical eddyviscosit ycoecien tis Preprin tsubmitted toElsevier Science 31July2003 z=D0,SeabedSeasurface z=0q(z) yx z(z)=A v(z)dq=dz(0)=w(cos()^x+sin()^y) A v(z),eddy viscosit y Fig.1.Thegeneral physical modelofthedepth-dep enden teddyviscosit yoceanog- raphyproblem represen tedasacontinuously dieren tiable function ofdepth. Ourgoalisto construct theblockmatrix thatrepresen tsthisSinc-Galerkin formulation and useitforsolving acoupled linear system. 2Governing Equations Todevelopamathematical model,werstconstruct aright-handed coordi- natesystem withthevertical coordinatezdirected positivedownwardfrom thefreesurface, andwithxandydirected northwardandeastward,respec- tively.Aplane atz=D0corresp ondstotheimpermeable boundary atthe seabed.Forthepurposeofillustrating theSinc-Galerkin metho d,weemplo y severalofthesimplifying assumptions invokedinone-dimensional wind-drift studies -theoceandepth,D0,andoceanmass densit y,,areassumed con- stant,andtheeects oftides, inertial terms, freesurface slope,andvariations inatmospheric pressure areneglected. Curren tsaredrivenbyatangen tial surface wind stress (thesurface isatz=0)ofmagnitude w,represen ted as(0)=w(cos()^x+sin()^y);withbeingtheangle betweenthepos- itivex-axisandthewind direction inFigure 1.Here ^xand^yrepresen t unitvectors inthedirections ofthepositivex-axisandthepositivey-axis, respectively.Inthismodel,called aspecied eddy viscosit ymodel,internal frictional stresses areparameterized as(z)=A v(z)dq=dz,where the specied eectiv evertical eddyviscosit ycoecien tA v(z)isacontinuously dieren tiable function ofz2(0;D0).Hereq(z)=U(z)^x+V(z)^y represen tsthehorizon talwind-drift curren twhichisthedierence between thetotalvelocityandthegeostrophic curren tin[2].TheCoriollis forceacting onthemovingwaterisbalanced byahorizon talpressure gradien tforce. 2 Under thepresen tassumptions, theconserv ation oflinear momen tumequa- tions express abalance betweentheCoriollis forceandtheinternal friction associated withturbulence. Thewind-drift curren tqisdetermined bysolving theboundary-v alueproblem d dz A v(z)dq dz! =f^zq;0<z<D0; (1) where thestress condition attheseasurface,z=0,isthetangen tialsurface windstress A v(0)dq(0) dz=w(cos()^x+sin()^y) (2) while attheseabed,z=D0,thefrictional stress islinearly proportional to thecurren t,hence A v(D0)dq(D0) dz=kfq(D0): (3) Heref2 sinistheCoriollis parameter atlatitudewhere =7:29 105rads1istheangular speedofrotation oftheearth. Though =2<< =2,fortherestofourworkweassume weareinthenorthern hemisphere and hence 0<<=2(f>0).Theparameterkfisthelinear slipbottom stress coecien t.Ifq(D0)=0,(3)iscalled azero-v elocity(no-slip) condition. From(1),wecansimplify theboundary-v alueproblem asfollows d dz A v(z)dq dz! =d dz A v(z)dU(z) dz! ^x+d dz A v(z)dV(z) dz! ^y =f^zq =f^z[U(z)^x+V(z)^y] =f U(z)^yV(z)^x ;0<z<D0: Then itcanbewritten incomponentformas d dz A v(z)dU(z) dz! =fV(z);0<z<D0 (4) and d dz A v(z)dV(z) dz! =fU(z);0<z<D0: (5) Thestress condition attheseasurface,z=0,separates as A v(0)dU(0) dz=wcos();A v(0)dV(0) dz=wsin(): (6) 3 Attheseabed,z=D0,thefrictional stress separates as A v(D0)dU(D0) dz=kfU(D0);A v(D0)dV(D0) dz=kfV(D0):(7) Tonondimensionalize themodelequations webeginwithameasure ofnear surface turbulen teddyviscosit y,A0A v(0),anddene anominal \upper- layer"Ekman depth byDEq 2A0=f.Alsodene acurren tspeedinunitsof U0=wDE=(A0)=p 2w=(pA0f).U0isthenatural velocityscaleinanin- nitely deepseawithuniform eddyviscosit yinthesteady-state. Anondimen- sional formoftheequations ofmotion canbeexpressed withtheintroduction ofthenondimensional variables zz D0;Av(z)A v(z) A v(0);q(z)q(z) U0U(z)^x+V(z)^y (8) together withtwonondimensional constan ts(adepth ratioandabottom friction parameter)givenby D0 DE=D0s f 2A0;A0Av(1) kfD0=A v(D0) kfD0: (9) Weapply (8)and(9)tonondimensionalize theequations (4)and(5).For(4), wehave 1 D0d dz A0Av(z)1 D0d(U0U(z)) dz! =fU0V(z);0<z<1; whichresults in d dz Av(z)dU(z) dz! =22V(z);0<z<1: (10) Similarly ,theresult from(5)is d dz Av(z)dV(z) dz! =22U(z);0<z<1: (11) Thesurface boundary conditions in(6)arenondimensionalized toproduce A0 D0d(U0U(0)) dz=wcos();A0 D0d(U0V(0)) dz=wsin(); whichleads to dU(0) dz=cos();dV(0) dz=sin(): (12) 4 Theseabedconditions in(7)arealsonondimensionalized andbecome U(1)+dU(1) dz=0;V(1)+dV(1) dz=0: (13) Inthespecialno-slip case,wecanset=0in(13). Wersttransform thenonhomogenous boundary conditions tohomogeneous boundary conditions byusing thelinear transformations U(z)=u(z)+(1+z)cos();V(z)=v(z)+(1+z)sin():(14) Therstderivativeofeachtransformation in(14)yields dU(z) dz=du(z) dzcos();dV(z) dz=dv(z) dzsin(); sotheresulting boundary-v alueproblem foru(z)satises d dz Av(z)du dz! +d dz(Av(z)cos())=22[v(z)+(1+z)sin()]; 0<z<1; whichgives d dz Av(z)du dz! +cos()A0 v(z)=22v(z)23(1+z)sin(); 0<z<1: (15) Similarly theboundary-v alueproblem forv(z)satises d dz Av(z)dv dz! +sin()A0 v(z)=22u(z)+23(1+z)cos(); 0<z<1: (16) Thesurface boundary conditions (12)become du(0) dz=0;dv(0) dz=0 (17) while theseabedboundary conditions (13)become u(1)+du(1) dz=0;v(1)+dv(1) dz=0: (18) 5 Forthepurposeofillustrating theexposition oftheSinc-Galerkin technique, wedene Lu(z)d dz Av(z)du dz! ;Lv(z)d dz Av(z)dv dz! : (19) Then (15)and(16)arenowgivenbythecoupleduandvequation systems Lu(z)+22v(z)=F1(z);0<z<1; (20) Lv(z)22u(z)=F2(z);0<z<1; (21) where F1(z)=23(1+z)sin()cos()A0 v(z); (22) F2(z)=23(1+z)cos()sin()A0 v(z): (23) Thesurface boundary conditions are du(0) dz=0;dv(0) dz=0: (24) Theseabedboundary conditions are u(1)+du(1) dz=0;v(1)+dv(1) dz=0: (25) 3TheSincFunction andtheSinc-Galerkin Metho d Werstreview sincfunction properties, sincquadrature rules, andtheSinc- Galerkin metho d.These arediscussed thoroughly in[3]and[4].Thesinc function isdened forallz2Cby sinc(z)8 >>>>>< >>>>>:sin(z) z;ifz6=0 : 1;ifz=0 Forh>0andk=0;1;2;3;:::thetranslated sincfunctions withevenly spaced nodesaregivenby 6 S(k;h)(z)sinc zkh h! 8 >>>>>>< >>>>>>:sin(zkh h) zkh h;ifz6=kh : 1;ifz=kh TheSinc-Galerkin procedure forthecoupled linear system in(20)-(21) begins byselecting compositesincfunctions appropriate totheinterval(0;1)asthe basis functions fortheexpansion ofapproximate solutions forthecurren t componentsu(z)andv(z).Weintroducetheconformal mapping function w=(z)=lnz 1z ; (26) whichisaconformal mapping fromDE,theeye-shap eddomain inthezplane, ontotheinnite stripinthew-plane, DS,where DE=( z=x+iy:argz 1z<d 2) ; DS= w=u+iv:jvj<d 2 : x zplane DEdiy 0 1iv wplane u DS 1d Fig.2.Therelationship betweentheeye-shap eddomain, DE,andtheinnite strip, DS ThisisshowninFigure 2.Thebasisfunctions arederivedfromthecomposite translated sincfunctions, Sk(z)S(k;h)(z)=sinc (z)kh h! ; (27) forz2DE.These areshowninFigure 3forrealvalues,x. 7 PSfrag replacemen ts 0 00.1 0.3 0.5 0.7 0.91 10.8 0.80.6 0.60.4 0.40.2 0.2-0.4-0.2 xS(k;h)(x)Sincbasison(0;1) k=1 k=0 k=1 Fig.3.Three adjacen tmembersS(k;h)(x)when k=1;0;1andh= 8ofthe mappedsincbasisontheinterval(0;1) FortheNeumann orradiation boundary conditions (24)-(25), thesincbasis functions in(27)donothaveaderivativewhenztends to0or1.Thuswe modifythesincbasisfunctions as S(k;h)(z) 0(z)sinc(z)kh h 0(z): (28) These areshowninFigure 4forrealvalues,x.Notethatthederivativesofthe modied sincbasisfunctions aredened aszapproac hes0or1.Wealsoadd boundary basisfunctions thatareHermite polynomials. Thusthederivatives atz=0and1forthesebasisfunctions aredened. These Hermite polynomials aregivenby B0(z)=(2z+1)(1z)2;B1(z)=(1z)z2+(32z)z2:(29) SincrulesforaspecialclassoffunctionsB(DE)havebeendeveloped.Adis- cussion oftheproperties offunctions inB(DE)isfound in[3]and[4]. Denition 1Let:DE!DSbeaconformal mapofDEtoDSwithinverse .Let=f (u)2DE:1<u<1g=(0;1).ThenB(DE)istheclassof functionsFwhich areanalytic inDE,satisfy Z (t+L)jF(z)jdz!0;t!1; whereL=n iv:jvj<d 2o ,andontheboundary ofDE,denote d@DE, satisfy N(F)Z @DEjF(z)jdz<1: 8 PSfrag replacemen ts 0 00.1 0.3 0.5 0.7 0.9 1 0.8 0.6 0.40.2 -0.1-0.050.050.1 0.10.150.2 0.20.25 x(S(k;h)(x))=0(x)Sincbasison(0;1) k=1 k=0 k=1 Fig.4.Three adjacen tmembersS(k;h)(x) 0(x)when k=1;0;1andh= 8of themodied sincbasisontheinterval(0;1) Theproofsofthefollowingtheorems arefound in[3]and[4]. Theorem 1If0F2B(DE)thenforallz2; F(z)=1X j=1F(zj)S(j;h)(z)+EF where EFsin (z) h! 2iZ @DE0(w)F(w)dw ((w)(z))sin((w)=h): Theorem 2IfF2B(DE)then Z1 0F(z)dz=h1X j=1F(zj) 0(zj)+IF where IFi 2Z @DEF(z)(;h)(z) sin((z)=h)dz with (;h)(z)=exp"i(z) hsgn(I((z)))# : Theinnite quadrature ruleappearing inTheorem 2canbeevaluated di- rectly,butingeneral itmustbetruncated toanitesumfortheSinc-Galerkin metho d.Thefollowingtheorem indicates theconditions under whichexponen- tialconvergence results. 9 Theorem 3LetF2B(DE)andbeaconformal mapwithconstants,, andCsothat F(z) 0(z)C8 >>>>>< >>>>>:exp(j(z)j);z2a exp(j(z)j);z2b where afz2:(z)=x2(1;0)g= 0;1 2! ; bfz2:(z)=x2[0;1)g="1 2;1! : Thenthesinctrapezoidal quadratureruleis Z1 0F(z)dz=hNX j=MF(zj) 0(zj)+O(exp(Mh))+O(exp(Nh)) +O(exp(2d=h)): (30) Hence,maketheselections N=" M+1# ;h=q 2d=(M); where[j
j]denotes thegreatestinteger,andtheexponential orderofthesinc trapezoidal quadraturerulein(30)isO(exp(p 2dM)). Corollary 1Animportant specialcaseofTheorem3occurswhentheinte- grandhastheformG(z)S(l;h)(z).Duetotheinterpolation S(l;h)(zj)=S(l;h)(jh)=(0) jl; thesincquadratureruleisaweighte dpointevaluation totheorderofthe method, Z1 0G(z)S(l;h)(z)dz=hG(zl) 0(zl)+O(exp(2d=h)): (31) IntheSinc-Galerkin technique, werstassume theapproximate solutions for u(z)andv(z)in(20)-(21), subjecttothemixed conditions (24)and(25),are represen tedby 10 ua(z)=cN1B0(z)+uh(z)+cN+1B1(z); (32) va(z)=dN1B0(z)+vh(z)+dN+1B1(z); (33) where choosingM=Ngives uh(z)=NX j=NcjSj(z) 0(z);vh(z)=NX j=NdjSj(z) 0(z): (34) Bothapproximate solutions clearly satisfy theboundary conditions (24)-(25). Them=2N+3coecien tsn cjoN+1 j=N1andthemcoecien tsn djoN+1 j=N1 aredetermined byorthogonalizing theresidual Lua(z)+22va(z)F1(z)and Lva(z)22ua(z)F2(z)withrespecttothesincbasisfunctionsn SjoN+1 j=N1 in(27).Theinner product (F;G)=Z1 0F(z)G(z)w(z)dz; usesaweightfunctionw(z)=1=q 0(z)=q z(1z).Sothisyields the discrete Sinc-Galerkin system, forj=N1;:::;N+1  Lua+22vaF1;Sj =0; (35)  Lva22uaF2;Sj =0: (36) Thisleads to cN1(LB0;Sj)+(Luh;Sj)+cN+1(LB1;Sj)+ 22va;Sj =(F1;Sj);(37) dN1(LB0;Sj)+(Lvh;Sj)+dN+1(LB1;Sj)+ 22ua;Sj =(F2;Sj):(38) With thesincquadrature rulein(31)weevaluate theinner products in(37)- (38).First, weevaluate (LB0;Sj)and(LB1;Sj)as (LB0;Sj)=Z1 0LB0(z)Sj(z)q 0(z)dzhLB0(zj) (0(zj))3=2; (39) (LB1;Sj)=Z1 0LB1(z)Sj(z)q 0(z)dzhLB1(zj) (0(zj))3=2: (40) Thenodalpointszj=ejh=(1+ejh).Thenumerators inthese terms are 11 LB0(z)=d dz Av(z)dB0(z) dz! =6d dz Av(z)z(1z)! =6 Av(z)(12z)+A0 v(z)z(1z)! ; LB1(z)=d dz Av(z)dB1(z) dz! =d dz Av(z)(23z+66z)z! = Av(z)(2+66z12z)+A0 v(z)(23z+66z)z! : Theinner products (Luh;Sj)and(Lvh;Sj)in(37)-(38) areevaluated byin- tegrating byparts twice. Thussince Luh(z)=NX k=Nck0 @Av(z) Sk(z) 0(z)!01 A0 ; (41) theinner productinvolving Luhis (Luh;Sj)=NX k=Nck0 B@0 @Av Sk 0!01 A0 ;Sj1 CA: (42) Theinner productisevaluated byintegrating byparts twicetoget 0 B@0 @Av Sk 0!01 A0 ;Sj1 CA=Z1 00 @Av(z) Sk(z) 0(z)!01 A0 Sj(z)q 0(z)dz (43) =BTZ1 0Sk(z) 0(z)Av(z)0 @Sj(z)q 0(z)1 A00 dz Z1 0Sk(z) 0(z)A0 v(z)0 @Sj(z)q 0(z)1 A0 dz; where BT=Av(z) Sk(z) 0(z)!0S(z)q 0(z)1 0+Av(z)Sk(z) 0(z)0 @Sj(z)q 0(z)1 A01 0: BTtends tozeroasshownin[1]. Since(z)=lnz 1z ,then1 0(z)=z(1z).Thisleads totheresult 12 1q 0(z) 1q 0(z)!0 =1 2(12z),1 (0(z))3=2 1q 0(z)!00 =1 4:Thederivative in(43)are 0 @Sj(z)q 0(z)1 A0 =dSj(z) dq 0(z)+Sj(z) 1q 0(z)!0 = dSj(z) d+1 2(12z)Sj(z)!q 0(z); (44) and 0 @Sj(z)q 0(z)1 A00 =d2Sj(z) d2(0(z))3=2+dSj(z) d 00(z)q 0(z)+20(z)00(z) 2 0(z)3=2! +Sj(z) 1q 0(z)!00 = d2Sj(z) d21 4Sj(z)! (0(z))3=2: (45) Applying thesincquadrature rule,(43)becomes 0 B@0 @Av Sk 0!01 A0 ;Sj1 CA=Z1 0Sk(z) 0(z)Av(z) d2Sj(z) d21 4Sj(z)! (0(z))3=2dz Z1 0Sk(z) 0(z)A0 v(z) dSj(z) d+1 2(12z)Sj(z)!q 0(z)dz h1 h2(2) jk+1 4(0) jk Av(zk)(0(zk))1=2 +h1 h(1) jk+1 2(2zk1)(0) jk A0 v(zk)(0(zk))3=2: So (Luh;Sj)NX k=N1 h2(2) jk+1 4(0) jk Av(zk)(0(zk))1=2ck (46) +NX k=N1 h(1) jk+1 2(2zk1)(0) jk A0 v(zk)(0(zk))3=2ck: HeretherthderivativeofS(j;h)(z),withrespectto,evaluated atthe 13 nodalpointzkisdenoted by 1 hr(r) jkdr dr[S(j;h)(z)] z=zk: (47) Similarly (Lvh;Sj)NX k=N1 h2(2) jk+1 4(0) jk Av(zk)(0(zk))1=2dk (48) +NX k=N1 h(1) jk+1 2(2zk1)(0) jk A0 v(zk)(0(zk))3=2dk: Next, applying thesincquadrature rule(31),theinner products (ua;Sj)and (va;Sj)are (ua;Sj)=Z1 0ua(z)Sj(z)q 0(z)dzh1 (0(zj))3=2ua(zj); (49) (va;Sj)=Z1 0va(z)Sj(z)q 0(z)dzh1 (0(zj))3=2va(zj); (50) where thesolutionsua(zj)andva(zj)aregivenby ua(zj)=cN1B0(zj)+cj1 0(zj)+cN+1B1(zj); (51) va(zj)=dN1B0(zj)+dj1 0(zj)+dN+1B1(zj): (52) Last, applying thesincquadrature rule(31),theinner products (F1;Sj)and (F2;Sj)are (F1;Sj)=Z1 0F1(z)Sj(z)q 0(z)dzh1 (0(zj))3=2F1(zj); (53) (F2;Sj)=Z1 0F2(z)Sj(z)q 0(z)dzh1 (0(zj))3=2F2(zj); (54) whereF1andF2aregivenin(22)-(23). Substituting (39),(40),(46),(48),(51), (52),(53),and(54)into(37)-(38), theresult willbelinear equations whose solutions givethecoecien tsforua(z)andva(z).ForN1jN+1; 14 LB0(zj) (0(zj))3=2cN1+NX k=N1 h2(2) jk+1 4(0) jk Av(zk)(0(zk))1=2ck +NX k=N1 h(1) jk+1 2(2zk1)(0) jk A0 v(zk)(0(zk))3=2ck (55) +LB1(zj) (0(zj))3=2cN+1+22(0(zj))3=2va(zj)=F1(zj) (0(zj))3=2 and LB0(zj) (0(zj))3=2dN1+NX k=N1 h2(2) jk+1 4(0) jk Av(zk)(0(zk))1=2dk +NX k=N1 h(1) jk+1 2(2zk1)(0) jk A0 v(zk)(0(zk))3=2dk (56) +LB1(zj) (0(zj))3=2dN+122(0(zj))3=2ua(zj)=F2(zj) (0(zj))3=2: Introducethecolumn vectors canddthatrepresen tthecoecien tsofthe approximate solutions in(32)-(33), c[cN1cN


cNcN+1]T;d[dN1dN


dNdN+1]T; anddene F1[F1(zN1)F1(zN)


F1(zN)F1(zN+1)]T; F2[F2(zN1)F2(zN)


F2(zN)F2(zN+1)]T; ua[ua(zN1)ua(zN)


ua(zN)ua(zN+1)]T; va[va(zN1)va(zN)


va(zN)va(zN+1)]T; whereF1(z)andF2(z)aregivenin(22)-(23). Theexpressions in(47)foreach jandkcanbestored inamatrixI(r)=[(r) jk].Forr=0;1;2,wehave (0) jk[S(j;h)(z)]jz=zk=8 >>>>>< >>>>>:1;ifk=j 0;ifk6=j; (57) 15 (1) jkhd d[S(j;h)(z)]jz=zk=8 >>>>>< >>>>>:0;ifk=j (1)kj kj;ifk6=j; (2) jkh2d2 d2[S(j;h)(z)]jz=zk=8 >>>>>< >>>>>:2 3;ifk=j 2(1)kj (kj)2;ifk6=j: Thefollowingmatrices willbesome examples forI(0);I(1);I(2).GivenN 1jN+1,(m=2N+3),themm,square matricesI(0);I(1);I(2)are givenby I(0)=2 6666641:::0 ......... 0:::13 777775=I;I(1)=2 6666666666664011 2:::(1)m1 m1 1............ 1 2.........1 2 ............1 (1)m m1:::1 2103 7777777777775; (58) I(2)=2 66666666666642 322 22:::2(1)m1 (m1)2 2............ 2 22.........2 22 ............2 2(1)m1 (m1)2:::2 2222 33 7777777777775: (59) When NjN,weremovetherstandlastcolumns ofI(0),I(1),and I(2)in(58)-(59), toarriveatthemn,(n=2N+1),non-square matrices I(0) z=2 66666666666640:::0 1:::0 ......... 0:::1 0:::03 7777777777775;I(1) z=2 666666666666666411 2:::(1)m2 m2 0......... 1......1 2 1 2......1 .........0 (1)m1 m2:::1 213 7777777777777775; (60) (61) 16 I(2) z=2 66666666666666642 2 22:::2(1)m2 (m2)2 2 3......... 2......2 22 2 22...... 2 .........2 3 2(1)m2 (m2)2:::2 22 23 7777777777777775: (62) Thennsquare diagonal matrix Dn(f)iswritten as Dn(f)=2 6666666666664f(zN) 0 ... f(z0) ... 0 f(zN)3 7777777777775: (63) From(55)and(56)wearriveatthediscrete system Bbc+22Dm 1 (0)3=2! Ebd=Dm 1 (0)3=2! F1 Bbd22Dm 1 (0)3=2! Ebc=Dm 1 (0)3=2! F2; where Bbisdened bythemmbordered matrix Bb[aN1jAnsjaN+1]: (64) Thenon-square mnmatrix Ans1 h2I(2) z+1 4I(0) z Dn Av (0)1=2! (65) +1 hI(1) z+1 2I(0) zDn(2z1) Dn A0 v (0)3=2! andthem1column vectors aN1andaN+1havejthcomponent, N1jN+1; [aN1]jLB0(zj) (0(zj))3=2;[aN+1]jLB1(zj) (0(zj))3=2: (66) 17 Themmevaluator matrix Ebisdened by Eb" bN1jI(0)Dn 1 0! jbN+1# (67) andthem1column vectors bN1andbN+1havejthcomponent [bN1]jB0(zj);[bN+1]jB1(zj): (68) Thisleads tothecoupled linear discrete system AX=C (69) where the(2m)(2m)blockmatrix (m=2N+3) A=2 666666664Bb 22Dm 1 (0)3=2 Eb 22Dm 1 (0)3=2 Eb Bb3 777777775(70) andthe(2m)1column vectors X=2 666666664c d3 777777775;C=2 666666664Dm 1 (0)3=2 F1 Dm 1 (0)3=2 F23 777777775: (71) Itfollowsthatthesolutions foruaandvaatthenodalpointsaretheele- mentsofthiscolumn vectorXthatarecomputed fromthe(2m)(2m)block evaluator matrix2 666666664ua va3 777777775=2 666666664EbO O Eb3 7777777752 666666664c d3 777777775; (72) where Oisanmmzeromatrix. These solutions arerelated totheap- proximate nondimensional curren tcomponentsevaluated atthenodalpoints as Ua(zj)=ua(zj)+(1+zj)cos(); Va(zj)=va(zj)+(1+zj)sin(): 18 4Numerical Testing: Constan tEddy Viscosit y Weillustrate theaccuracy oftheSinc-Galerkin metho dwhen applied tothe constan teddyviscosit ymodelformulated asasystem ofreal-valued coupled ordinary dieren tialequations. Since thegoverning equations andvariables werenondimensionalized, theonlyoperativeconstan tsin(20)-(25) are;, and.Inrelating these parameters tothe\constan tsofnature", weadopt the followingnominal values:f=0:0001s1(appropriate totemperate northern latitudes), seawaterdensit y=1103kgm3,andairdensit yair= 1:25kgm3. Surface windstress isassumed toberelated tothesquare ofthewindspeed Ww(inms1)by w=CDairW2 w; (73) where thedimensionless parameterCD0:0012forWw<12ms1,there- afterincreasing linearly toabout0.0025 atgaleforcewinds (Ww30ms1) [5].Inkeeping with[6]and[7],thelinear slipbottom stress coecien tkfis assigned avalueof0:002ms1forcomparison withother workandtodrama- tizethechange incurren tspeedoverthewatercolumn. Inpractice, avalueof kflowerbyanorder ofmagnitude maybepreferred. Field evidence suggests thatthenear-surface valueofthevertical eddyviscosit yisrelated tothewind speed.Carter [8]hassuggested that,ifthewindisnotfetch-limited andthe seastateisfullydeveloped,thenA v(0)inunits ofm2s1isgivenby A v(0)0:304104W3 w: (74) With theparameters andrelationships above,andkeeping inmind ourdesire tocompare results withthose in[9],wechooseourconstan teddyviscosit yto be A v(z)0:02m2s1(75) withw=p 2=10=0:1414Nm2.SinceDE=q 2Av(0)=f=20m,wethen have=D0=DE=D0=20.Inkeeping with[9],wewilluseD0=100mand hence=5,whichitwillbethroughout. Thenumerical results arecompared totheexact solutionW(z)=U0[U(z)+ iV(z)]whereU(z)isgivenby U(z)=R(Wc(z))cos()I(Wc(z))sin() (76) andV(z)isgivenby V(z)=R(Wc(z))sin()+I(Wc(z))cos(): (77) HereR(Wc(z))andI(Wc(z))denote therealandimaginary parts ofWc(z), 19 respectively,and Wc(z)=(1i)cosh((1i)(1z))+sinh((1i)(1z)) (1i)[cosh ((1i))+(1i)sinh((1i))]: Theresults oftheSinc-Galerkin approximationsUa(zj)andVa(zj)werecom- pared withtheexact solutions forU(zj)andV(zj)atthesincgridpointsS withh==p 2Ngivenby S=n zj=1(jh)=ejh=(ejh+1):j=N1;:::;N+1o : (78) These results werethenmultiplied bythenatural velocityscaleU0togivea dimensional represen tation ofthevelocities. Allnumerical simulations wererunonaSUNBLADE 1000withMATLAB Version 6.1.Toillustrate theperformance ofthemetho d,themaxim umabso- luteerrors arereported as kUSk= max N1jN+1fU0jUa(zj)U(zj)jg; kVSk= max N1jN+1fU0jVa(zj)V(zj)jg; and kESk=maxfkUSk;kVSkg; (79) where theunits arems1. Throughout, comparable graphs (ofeddyviscosit yfunctions andvelocitycom- ponents)areshownonthesamescale. Thehorizon talprojections oftheEkman spirals arealsoshownonthesamescale. Thiswayvisual comparisons ofthese various quantities arereadily made. Example 1(Linear stress condition attheseabed)Forthisexample we choose=45oandforthelinear stress condition attheseabedwehave =A v(D0)=(kfD0)=0:1.Wendtheapproximate solutionUa(z)andVa(z), respectively,using thecoupled discrete system ofsize2m2m(m=2N+3), givenin(69)whereU0Ua(z)=U a(z)andU0Va(z)=V a(z).Theerrors are giveninTable1andillustrate theclassic exponentialconvergence typical of Sinc-Galerkin metho ds. Figure 5graphically depicts thenumerical convergence oftheSinc-Galerkin metho dtotheexact solution withthelinear stress bottom boundary condition, asNisrepeatedly doubled insize.Thehorizon talprojection oftheEkman spiral forN=64isindistinguishable fromthetruesolution shownasthesolid line. 20 N2m hkUSk(ms1)kVSk(ms1)kESk(ms1) 4221.111 1.080e-03 7.478e-04 1.080e-03 8380.785 2.496e-04 1.293e-04 2.496e-04 16 700.555 2.762e-05 1.302e-05 2.762e-05 32134 0.393 8.985e-07 4.240e-07 8.985e-07 64262 0.278 5.775e-09 2.770e-09 5.775e-09 128 518 0.196 4.053e-12 1.937e-12 4.053e-12 Table1 Errors forExample 1(constan teddyviscosit y)onthesincgridSwiththelinear stress bottom condition for=0:1;=45o;=5;D0=100m;DE=20m −0.005PSfrag replacemen ts -0.02 0 0 0.02 0.04 0.06 0.08 0.1 0.12-0.04-0.035-0.03-0.025-0.02 -0.02-0.015-0.010.0050.01N=4 N=8 N=16 N=32 N=64 TrueNorth wardcurren tcomponentU a(m/s) Eastwardcurren tcomponentV a(m/s)Sinc-Galerkin Ekman spiral projections forincreasingN Fig.5.Sinc-Galerkin Ekman spiral projections forExample 1withincreasing Nforconstan teddyviscosit ywiththelinear stress bottom boundary condition for =0:1;=45o;=5;D0=100m;DE=20m Example 2(No-slip condition attheseabed)Weset=0and=45o andsolveforapproximate solutionsUa(z)andVa(z),respectively,byusing thecoupled discrete system in(69).Thereported errors areshowninTable2 andareverysimilar tothose forExample 1.Theaccuracy isnodieren tfor theNeumann condition inthisexample thanitisforthemixed condition in Example 1.Thehorizon talprojection oftheEkman spirals isshowninFigure 6.Thisgure portraystheexponentialconvergence ofthemetho d. 21 N2m hkUSk(ms1)kVSk(ms1)kESk(ms1) 4221.111 1.071e-03 7.500e-04 1.071e-03 8380.785 2.483e-04 1.297e-04 2.483e-04 16 700.555 2.752e-05 1.306e-05 2.752e-05 32134 0.393 8.957e-07 4.255e-07 8.957e-07 64262 0.278 5.756e-09 2.781e-09 5.756e-09 128 518 0.196 4.073e-12 1.975e-12 4.073e-12 Table2 Errors forExample 2(constan teddyviscosit y)onthesincgridSwiththezero- velocitybottom condition for=0;=45o;=5;D0=100m;DE=20m −0.005PSfrag replacemen ts -0.02 0 0 0.02 0.04 0.06 0.08 0.1 0.12-0.04-0.035-0.03-0.025-0.02 -0.02-0.015-0.010.0050.01N=4 N=8 N=16 N=32 N=64 TrueNorth wardcurren tcomponentU a(m/s) Eastwardcurren tcomponentV a(m/s)Sinc-Galerkin Ekman spiral projections forincreasingN Fig.6.Sinc-Galerkin Ekman spiral projections forExample 2withincreasing N forthecaseofconstan teddyviscosit ywithno-slip bottom boundary condition for =0;=45o;=5;D0=100m;DE=20m 5Numerical Testing: Variable Eddy Viscosit y Thefollowingexamples showtheresults forthesystem ofreal-valued cou- pledordinary dieren tialequations witheither adecreasing orquadratic eddy viscosit yfunction givenby(80)or(81).Aneddy viscosit ywhichdecreases quadratically fromthevalueofA v(0)=0:02m2s1totheminim umvalueof 22 A v(D0)=0:00125 m2s1isgivenby A v(z)=0:02[1(:0075)z]2;0<z<D0=100m: (80) Agraph ofthisA v(z)isshowninFigure 7where itiscontrasted withthe constan teddyviscosit yofA v(z)0:02m2s1.Aquadratic eddyviscosit yPSfrag replacemen ts 0 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.110 20 30 40 50 60 70 80 90 100Decreasing A v(z) Constan tA v(z) Eddy viscosit yfunctionA v(z)m2s1Depthz(m)Decreasing eddyviscosit yfunction Fig.7.Eddy viscosit yfunctions A v(z)=0:02 1(:0075)z2m2s1and A v(z)0:02m2s1 function thatranges fromA v(0)=0:02m2s1backtoA v(D0)=0:02m2s1is givenby A v(z)=0:02[1+(:12)z(1(:01)z)];0<z<D0=100m;(81) whichisshowninFigure 8.Themodelwiththedecreasing andquadratic eddy viscosit yfunctions aresolvedbythecoupled discrete system in(69)ofsize 2m2m(m=2N+3).Alsotheresults fromthecoupled system scheme are compared totheresults in[1]whichusedacomplex-v alued ordinary dieren- tialequation. Example 3Forthisexample inaseaofdepthD0=100m,theparameters arechosen tobe=0:1,=45o,and(sinceDE=20m)=5.Thedecreas- ingeddyviscosit yfunctionA v(z)isgivenby(80).Wendtheapproximate solutionsUa(z)andVa(z)byusing thecoupled discrete system in(69).The results comparing theconstan teddy viscosit yandthedecreasing eddy vis- cosityareshowninFigure 9andFigure 10forN=32.Figure 9contrasts theEkman spiral projections while Figure 10depicts thevelocitycomponents U a(z)andV a(z)fromthesurface totheseabottom. Theconvergence of 23 PSfrag replacemen ts 0 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.110 20 30 40 50 60 70 80 90 100Quadratic A v(z) Constan tA v(z) Eddy viscosit yfunctionA v(z)(m2=s)Depthz(m)Quadratic eddyviscosit yfunction Fig.8.Eddy viscosit yfunctions A v(z)=0:02[1+(:12)z(1(:01)z)]m2s1and A v(z)0:02m2s1 theSinc-Galerkin Ekman spiral projections forthedecreasing eddyviscosit y function in(81)isillustrated inFigure 11.PSfrag replacemen ts -0.02 0 00.02 0.04 0.06 0.08 0.1 0.12-0.04-0.03-0.02 -0.02-0.010.010.02 0.02Ekman spiral projection forquadratic A v(z) Ekman spiral projection forconstan tA v(z)North wardcurren tcomponentU a(m/s) Eastwardcurren tcomponentV a(m/s)Sinc-Galerkin Ekman spiral projections Fig.9.Sinc-Galerkin Ekman spiral projections forExample 3withbothconstan t anddecreasing eddyviscosit yfunctions andlinear stress bottom boundary condition for=0:1;=45o;=5;D0=100m;DE=20m 24 PSfrag replacemen ts -0.1 -0.08 -0.06 -0.04 -0.020 0 0.02 0.04 0.06 0.08 0.110 20 30 40 50 60 70 80 90 100U afordecreasing A v(z) V afordecreasing A v(z) Uforconstan tA v(z) V aforconstan tA v(z) VelocitycomponentsU aandV a(m/s)Depthz(m)Sinc-Galerkin velocitycomponents Fig.10.Sinc-Galerkin northwardandeastwardcalculated velocityproles forEx- ample 3withconstan tanddecreasing eddy viscosit yfunctions andlinear stress bottom boundary condition for=0:1;=45o;=5;D0=100m;DE=20m −0.005PSfrag replacemen ts -0.02 0 0 0.02 0.04 0.06 0.08 0.1 0.12-0.04-0.035-0.03-0.025-0.02 -0.02-0.015-0.010.0050.01N=4 N=8 N=16 N=32 N=64 TrueNorth wardcurren tcomponentU a(m/s) Eastwardcurren tcomponentV a(m/s)Sinc-Galerkin Ekman spiral projections forincreasingN Fig.11.Sinc-Galerkin Ekman spiral projections forExample 3withincreasing Nfor thecaseofthedecreasing eddyviscosit yfunction andlinear stress bottom boundary condition for=0:1;=45o;=5;D0=100m;DE=20m Example 4With thequadratic eddyviscosit yfunction in(81),theparame- tersareagain chosen tobe=0:1;=45o,and=5.Theapproximate solutionsUa(z)andVa(z)areillustrated bygraphs comparing theconstan t eddyviscosit ywiththequadratic eddyviscosit yinFigure 12andFigure 13 forN=32.Theconvergence oftheSinc-Galerkin Ekman spiral projections 25 forthequadratic eddyviscosit yfunction in(81)isillustrated inFigure 14. Theapproximations arecalculated fromthe2m2mcoupled discrete system in(69).PSfrag replacemen ts -0.02 0 00.02 0.04 0.06 0.08 0.1 0.12-0.04-0.03-0.02 -0.02-0.010.010.02 0.02Ekman spiral projection forquadratic A v(z) Ekman spiral projection forconstan tA v(z)North wardcurren tcomponentU a(m/s) Eastwardcurren tcomponentV a(m/s)Sinc-Galerkin Ekman spiral projections Fig.12.Sinc-Galerkin Ekman spiral projections forExample 4withbothconstan t andquadratic eddyviscosit yfunctions andlinear stress bottom boundary condition for=0:1;=45o;=5;D0=100m;DE=20mPSfrag replacemen ts -0.1 -0.08 -0.06 -0.04 -0.020 0 0.02 0.04 0.06 0.08 0.110 20 30 40 50 60 70 80 90 100U aforquadratic A v(z) V aforquadratic A v(z) U aforconstan tA v(z) V aforconstan tA v(z) VelocitycomponentsU aandV a(m/s)Depthz(m)Sinc-Galerkin velocitycomponents Fig.13.Sinc-Galerkin northwardandeastwardcalculated velocityproles forExam- ple4withconstan tandquadratic eddyviscosit yfunctions andlinear stress bottom boundary condition for=0:1;=45o;=5;D0=100m;DE=20m 26 −0.005PSfrag replacemen ts -0.02 0 0 0.02 0.04 0.06 0.08 0.1 0.12-0.04-0.035-0.03-0.025-0.02 -0.02-0.015-0.010.0050.01 N=4 N=8 N=16 N=32 N=64 TrueNorth wardcurren tcomponentU a(m/s) Eastwardcurren tcomponentV a(m/s)Sinc-Galerkin Ekman spiral projections forincreasingN Fig.14.Sinc-Galerkin Ekman spiral projections forExample 4withincreasing Nfor thecaseofthequadratic eddyviscosit yfunction andlinear stress bottom boundary condition for=0:1;=45o;=5;D0=100m;DE=20m Next, wewanttocompare theapproximate solutionsUa(z)andVa(z)forthe complex velocitysystem in[1]andthenewcoupled system introduced here. Toillustrate thecomparativ eperformance ofbothmetho ds,themaxim umab- solute errors betweentheirrespectivenumerical approximations arereported as kUCk= max N1jN+1fU0jU1(zj)U2(zj)jg; kVCk= max N1jN+1fU0jV1(zj)V2(zj)jg; kECk=maxfkUCk;kVCkg; whereU1(z)andV1(z)arecomputed bythecomplex velocitysystem in[1] andU2(z)andV2(z)arecomputed bythecoupled system in(69). Example 5Thisexample compares thesolution introduced heretothatcom- puted bythecomplex velocitysystem in[1].Theexample usedisthatde- scribedinExample 1.Allcomparisons showasimilarit ygreater thanthe accuracy ofthemetho dreported inTable1 27 N m2m kUCk kVCk kECk 411 221.005e-15 8.163e-16 1.005e-15 819 389.419e-16 1.193e-15 1.193e-15 16 35 706.279e-16 6.217e-15 6.217e-15 32 67134 1.130e-15 2.009e-15 2.009e-15 64131 262 5.651e-15 1.155e-14 1.155e-14 128 259 518 1.463e-14 1.758e-14 1.758e-14 Table3 Comparison betweentheapproximate solutions forthecomplex velocitysystem andthecoupled system byusing thesame sincgridsizeforExample 1forthe caseofconstan teddy viscosit ywithlinear stress bottom boundary condition for =0:1;=45o;=5;D0=100m;DE=20m Example 6Thisexample compares thesolution introduced heretothecom- plexvelocitysystem in[1].Thecomparison usedExample 4.Again theresults areextremely similar. Notethatthemetho dusedhererequires noneofthe complex arithmetic necessary intheworkof[1]. N m2m kUCk kVCk kECk 411 221.884e-16 9.419e-16 9.419e-16 819 383.768e-16 1.193e-15 1.193e-15 16 35 702.198e-16 3.579e-15 3.579e-15 32 67134 1.381e-15 2.261e-15 2.261e-15 64131 262 4.969e-15 1.068e-14 1.068e-14 128 259 518 4.458e-14 2.003e-14 4.458e-14 Table4 Comparison betweentheapproximate solutions forthecomplex velocitysystem and thecoupled system onthesame sincgridforExample 4forthedecreasing eddy viscosit yA v(z)=:02(1+:12z(1:01z))withlinear stress bottom boundary condition for=0:1;=45o;=5;D0=100m;DE=20m 28 References [1]D.F.Winter,J.Lund, andK.L.Bowers.Wind-driv encurren tsinaseawitha variable eddyviscosit ycalculated byaSincfunction Galerkin technique. Internat. J.Numer. Metho dsFluids ,33:1041-1073, 2000. [2]J.Brown,A.Colling, D.Park,J.Phillips, D.Rothery ,andJ.Wright.Ocean Circulation, TheOpenUniversity,Keynes, 1989. [3]J.Lund andK.L.Bowers.SincMetho dsforQuadrature andDieren tial Equations, SIAM: Philadelphia, 1992. [4]F.Stenger. Numerical Metho dsBased onSincandAnalytic Functions, Springer- Verlag, NewYork,1993. [5]W.G.Large andS.Pond.Openocean momen tum uxmeasuremen tsin moderate tostrong winds. J.Phys.Oceanogr ,11:324-336, 1981. [6]A.M.DaviesandA.Owen.Three-dimensional numerical seamodelusing Galerkin metho dwithapolynomial basisset.Appl. Math. Modeling ,3:421-428, 1979. [7]N.S.Heaps. Onthenumerical solution ofthethree-dimensional hydrodynamical equations fortidesandstorm surges. Mem. Soc.Sci.Liege., Ser.6,1:143-180, 1971. [8]D.J.T.Carter. Estimation ofwavespectra fromwaveheightandperiod.Inst. Oceanogr. Sci.Rep,135,1982. [9]C.E.Naimie. ATurbulen tBoundary LayerModelfortheLinearized Shallo w WaterEquations, NUBBLE USER'S MANUAL (Release 1.1).Technical Report NML-96-1, Dartmouth College, July31,1996. 29