{VERSION 5 0 "IBM INTEL NT" "5.0" } {USTYLETAB {CSTYLE "Maple Input" -1 0 "Courier" 0 1 255 0 0 1 0 1 0 0 1 0 0 0 0 1 }{CSTYLE "2D Math" -1 2 "Times" 0 1 0 0 0 0 0 0 2 0 0 0 0 0 0 1 }{CSTYLE "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" -1 256 "" 1 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 1 14 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 264 "" 1 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 266 "" 1 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 267 "" 1 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Heading 1" -1 3 1 {CSTYLE "" -1 -1 "Times " 1 18 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 4 1 0 1 0 2 2 0 1 } {PSTYLE "Title" 0 18 1 {CSTYLE "" -1 -1 "" 1 18 0 0 0 0 0 1 1 0 0 0 0 0 0 0 }3 0 0 -1 12 12 0 0 0 0 0 0 19 0 }{PSTYLE "Author" 0 19 1 {CSTYLE "" -1 -1 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }3 0 0 -1 8 8 0 0 0 0 0 0 -1 0 }{PSTYLE "Normal" -1 256 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 2 2 2 2 2 2 1 1 1 1 }3 1 0 0 0 0 1 0 1 0 2 2 0 1 }} {SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT -1 40 "Partial Differential Equa tions PowerTool" }}{PARA 19 "" 0 "" {TEXT -1 16 "by Dr. Jim Herod" }}} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 257 69 "Section \+ 8.3 Numerical Solutions for the One Dimensional Wave Equation" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 0 {PARA 3 "" 0 "" {TEXT 264 34 " Maple Packages used in Section 8.3" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "with(Linea rAlgebra):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 27 "The problem we co nsider is " }}{PARA 0 "" 0 "" {TEXT -1 8 " " }{XPPEDIT 18 0 "di ff(u,`$`(x,2));" "6#-%%diffG6$%\"uG-%\"$G6$%\"xG\"\"#" }{TEXT -1 3 " = " }{XPPEDIT 18 0 "diff(u,`$`(t,2));" "6#-%%diffG6$%\"uG-%\"$G6$%\"tG \"\"#" }{TEXT -1 9 " - F(t,x)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 51 "with boundary conditions u(t, 0) = u (t,1) = 0," }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 52 "and with initial conditions u(0, x) = f(x) and " } {XPPEDIT 18 0 "diff(u,t);" "6#-%%diffG6$%\"uG%\"tG" }{TEXT -1 13 "(0,x ) = g(x)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 50 "By now, we expect to create a sequence of vectors " }{XPPEDIT 18 0 "U[1],U[2];" "6$&%\"UG6#\"\"\"&F$6#\"\"#" }{TEXT -1 6 ", ... " } {XPPEDIT 18 0 "U[M];" "6#&%\"UG6#%\"MG" }{TEXT -1 27 " along the time \+ axis. Each " }{XPPEDIT 18 0 "U[m];" "6#&%\"UG6#%\"mG" }{TEXT -1 57 " h as its values created from what has come before. These " }{TEXT 258 2 "U " }{TEXT -1 143 "'s are found from the preceding ones in a manner b ased on equations found by replacing the partial derivatives with appr opriate approximations." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 14 "If the vector " }{XPPEDIT 18 0 "U[m];" "6#&%\"UG6#% \"mG" }{TEXT -1 5 " has " }{TEXT 259 1 "N" }{TEXT -1 101 " terms, then the first and last term are determined by the boundary conditions. Th at is, for each m, " }{XPPEDIT 18 0 "U[m];" "6#&%\"UG6#%\"mG" }{TEXT -1 8 "(0) and " }{XPPEDIT 18 0 "U[m];" "6#&%\"UG6#%\"mG" }{TEXT -1 78 "(N) are zero for this problem. In a similar manner, we expect that th e vector " }{XPPEDIT 18 0 "U[0];" "6#&%\"UG6#\"\"!" }{TEXT -1 19 ", co rresponding to " }{TEXT 260 1 "t" }{TEXT -1 89 " = 0, is determined by the initial condition, u(0,x) = f(x). The implication is that, if " } {XPPEDIT 18 0 "x[n];" "6#&%\"xG6#%\"nG" }{TEXT -1 23 " is the mesh alo ng the " }{TEXT 261 1 "x" }{TEXT -1 11 " axes, then" }}{PARA 0 "" 0 " " {TEXT -1 11 " " }{XPPEDIT 18 0 "U[0];" "6#&%\"UG6#\"\"!" } {TEXT -1 8 "(n) = f(" }{XPPEDIT 18 0 "x[n];" "6#&%\"xG6#%\"nG" }{TEXT -1 24 ") for n = 1, 2, ... , N." }}{PARA 0 "" 0 "" {TEXT -1 94 "How th e remaining initial condition enters the approximation scheme will be \+ developed shortly." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 65 "We state the discrete replacement terms: Let h = 1/N and \+ k = 1/M." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 4 " " }{XPPEDIT 18 0 "diff(u,`$`(x,2));" "6#-%%diffG6$%\"uG-%\"$G6$ %\"xG\"\"#" }{TEXT -1 21 " is replaced by " }{XPPEDIT 18 0 "(U[m] (n+1)-2*U[m](n)+U[m](n-1))/(h^2);" "6#*&,(-&%\"UG6#%\"mG6#,&%\"nG\"\" \"F-F-F-*&\"\"#F--&F'6#F)6#F,F-!\"\"-&F'6#F)6#,&F,F-F-F4F-F-*$%\"hGF/F 4" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 3 "and" }}{PARA 0 "" 0 " " {TEXT -1 4 " " }{XPPEDIT 18 0 "diff(u,`$`(t,2));" "6#-%%diffG6$% \"uG-%\"$G6$%\"tG\"\"#" }{TEXT -1 28 " - F(t,x) is replaced by " } {XPPEDIT 18 0 "(U[m+1](n)-2*U[m](n)+U[m-1](n))/(k^2);" "6#*&,(-&%\"UG6 #,&%\"mG\"\"\"F+F+6#%\"nGF+*&\"\"#F+-&F'6#F*6#F-F+!\"\"-&F'6#,&F*F+F+F 46#F-F+F+*$%\"kGF/F4" }{TEXT -1 4 " - " }{XPPEDIT 18 0 "F[m](n);" "6# -&%\"FG6#%\"mG6#%\"nG" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 62 "From the statement of the problem, thes e should be equal. Let " }{XPPEDIT 18 0 "delta;" "6#%&deltaG" }{TEXT -1 3 " = " }{XPPEDIT 18 0 "h/k;" "6#*&%\"hG\"\"\"%\"kG!\"\"" }{TEXT -1 44 " . We let Maple solve this equality for us." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 92 "eq:=(U[m](n+1)-2*U[m](n)+U[m](n-1))/h^2 =\n (U[m+1](n)-2*U[m](n)+U[m-1](n))/k^2-F[m](n);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "solve(eq,U[m+1](n));" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 13 " Here is the " }{TEXT 262 17 "Summary Equation:" } }{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 2 " " } {XPPEDIT 18 0 "U[m+1](n);" "6#-&%\"UG6#,&%\"mG\"\"\"F)F)6#%\"nG" } {TEXT -1 6 " = " }{XPPEDIT 18 0 "delta^2*U[m](n-1)+2*(1-delta^2)*U[ m](n)+delta^2*U[m](n+1)-U[m-1](n)+k^2*F[m](n);" "6#,,*&%&deltaG\"\"#-& %\"UG6#%\"mG6#,&%\"nG\"\"\"F/!\"\"F/F/*(F&F/,&F/F/*$F%F&F0F/-&F)6#F+6# F.F/F/*&F%F&-&F)6#F+6#,&F.F/F/F/F/F/-&F)6#,&F+F/F/F06#F.F0*&%\"kGF&-&% \"FG6#F+6#F.F/F/" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 101 "This equation is appropriate for m > 1 a nd for 0 < n < N. We have already discussed what happens for " } {XPPEDIT 18 0 "U[0];" "6#&%\"UG6#\"\"!" }{TEXT -1 23 ". Here is the st ory on " }{XPPEDIT 18 0 "U[1];" "6#&%\"UG6#\"\"\"" }{TEXT -1 96 ". The re are two considerations. From the Summary Equation we see that what \+ should happen is that" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 3 " " }{XPPEDIT 18 0 "U[1](n);" "6#-&%\"UG6#\"\"\"6#%\"n G" }{TEXT -1 6 " = " }{XPPEDIT 18 0 "delta^2*U[0](n-1)+2*(1-delta^2 )*U[0](n)+delta^2*U[0](n+1)-U[-1](n)+k^2*F[0](n);" "6#,,*&%&deltaG\"\" #-&%\"UG6#\"\"!6#,&%\"nG\"\"\"F/!\"\"F/F/*(F&F/,&F/F/*$F%F&F0F/-&F)6#F +6#F.F/F/*&F%F&-&F)6#F+6#,&F.F/F/F/F/F/-&F)6#,$F/F06#F.F0*&%\"kGF&-&% \"FG6#F+6#F.F/F/" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 28 "What must be done with that " }{XPPEDIT 18 0 "U[-1];" "6#&%\"UG6#,$\"\"\"!\"\"" }{TEXT -1 89 " term comes from the remaining initial condition. The discrete replacement term we use is" }}{PARA 0 "" 0 "" {TEXT -1 11 " " }{XPPEDIT 18 0 "(U[1] (n)-U[-1](n))/(2*k);" "6#*&,&-&%\"UG6#\"\"\"6#%\"nGF)-&F'6#,$F)!\"\"6# F+F0F)*&\"\"#F)%\"kGF)F0" }{TEXT -1 14 " = g(n / N)." }}{PARA 0 "" 0 "" {TEXT -1 76 "We now know the replacement for the minus one term i n the equation defining " }{XPPEDIT 18 0 "U[1];" "6#&%\"UG6#\"\"\"" } {TEXT -1 1 ":" }}{PARA 0 "" 0 "" {TEXT -1 4 " " }{XPPEDIT 18 0 "U[1 ](n);" "6#-&%\"UG6#\"\"\"6#%\"nG" }{TEXT -1 6 " = " }{XPPEDIT 18 0 "delta^2*U[0](n-1)+2*(1-delta^2)*U[0](n)+delta^2*U[0](n+1)-(U[1](n)-2* k*g(n/N))+k^2*F[0](n);" "6#,,*&%&deltaG\"\"#-&%\"UG6#\"\"!6#,&%\"nG\" \"\"F/!\"\"F/F/*(F&F/,&F/F/*$F%F&F0F/-&F)6#F+6#F.F/F/*&F%F&-&F)6#F+6#, &F.F/F/F/F/F/,&-&F)6#F/6#F.F/*(F&F/%\"kGF/-%\"gG6#*&F.F/%\"NGF0F/F0F0* &FDF&-&%\"FG6#F+6#F.F/F/" }{TEXT -1 2 " ." }}{PARA 0 "" 0 "" {TEXT -1 3 "Or," }}{PARA 0 "" 0 "" {TEXT -1 4 " " }{XPPEDIT 18 0 "U[1](n);" "6#-&%\"UG6#\"\"\"6#%\"nG" }{TEXT -1 3 " = " }{XPPEDIT 18 0 "delta^2/2 *U[0](n-1)+(1-delta^2)*U[0](n)+delta^2/2*U[0](n+1)+k*g(n/N)+k^2*F[0](n )/2;" "6#,,*(%&deltaG\"\"#F&!\"\"-&%\"UG6#\"\"!6#,&%\"nG\"\"\"F0F'F0F0 *&,&F0F0*$F%F&F'F0-&F*6#F,6#F/F0F0*(F%F&F&F'-&F*6#F,6#,&F/F0F0F0F0F0*& %\"kGF0-%\"gG6#*&F/F0%\"NGF'F0F0*(F?F&-&%\"FG6#F,6#F/F0F&F'F0" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 116 "We write these replacement equations in the matrix form. For expo sition, suppose that N = 4 = M. Here is the matrix " }{TEXT 263 1 "A" }{TEXT -1 36 " that will be used in this analysis." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "N:=4; M:=4; \ndelta:=(1/N^2)/(1/M^2);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 145 "A:=Matrix([[0,0,0,0,0],[delta^2,1-delta^2,delta^2,0,0],\n [0, delta^2,1-delta^2,delta^2,0],\n [0,0,delta^2,1-delta^2,delta^2],[0 ,0,0,0,0]]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 11 "The vector " } {XPPEDIT 18 0 "U[0];" "6#&%\"UG6#\"\"!" }{TEXT -1 53 " is defined usin g only one of the initial conditions:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "U[0]:=Vector([f(0),f(1/4),f(2/4),f(3/4),f(1)]);" }}} {PARA 0 "" 0 "" {TEXT -1 11 "The vector " }{XPPEDIT 18 0 "U[1];" "6#&% \"UG6#\"\"\"" }{TEXT -1 124 " is defined using the above analysis. We \+ first take care of the non-homogeneous term, F(t,x), with a sequence \+ of vectors V." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 66 "for m from \+ 1 to M do\n V[m]:=Vector([seq(F(m/M,n/N),n=0..N)]):\nod:" }}}{PARA 0 "" 0 "" {TEXT -1 8 "Here is " }{XPPEDIT 18 0 "U[1];" "6#&%\"UG6#\"\"\" " }{TEXT -1 1 "." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "U[1]:=A. U[0]+k^2*V[1]+2*k*Vector([seq(g(n/N),n=0..N)]);" }}}{PARA 0 "" 0 "" {TEXT -1 14 "The remaining " }{TEXT 265 2 "U " }{TEXT -1 34 "'s are de termined by the iteration" }}{PARA 0 "" 0 "" {TEXT -1 16 " \+ " }{XPPEDIT 18 0 "U[m+1];" "6#&%\"UG6#,&%\"mG\"\"\"F(F(" }{TEXT -1 3 " = " }{XPPEDIT 18 0 "A*U[m]+V[m+1];" "6#,&*&%\"AG\"\"\"&%\"UG6#% \"mGF&F&&%\"VG6#,&F*F&F&F&F&" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 40 "A word should be said about the size of " }{XPPEDIT 18 0 "delta;" "6#%&deltaG" }{TEXT -1 188 " for wh ich the discrete replacement scheme can be expected to be appropriate. Reference is made to the Courant-Friedrichs-Lewy Convergence Criterio n. In this problem, the criteria is that " }}{PARA 0 "" 0 "" {TEXT -1 37 " 0 < " }{XPPEDIT 18 0 "delta <= 1; " "6#1%&deltaG\"\"\"" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 0 {PARA 3 "" 0 "" {TEXT 256 9 "Problem 1" }}{PARA 0 "" 0 "" {TEXT -1 29 "We work the following problem" }}{PARA 0 "" 0 "" {TEXT -1 8 " " }{XPPEDIT 18 0 "diff(u,`$`(x,2));" "6#-%%diffG6$%\"uG- %\"$G6$%\"xG\"\"#" }{TEXT -1 3 " = " }{XPPEDIT 18 0 "diff(u,`$`(t,2)); " "6#-%%diffG6$%\"uG-%\"$G6$%\"tG\"\"#" }{TEXT -1 10 " - 16 sin(" } {XPPEDIT 18 0 "Pi;" "6#%#PiG" }{TEXT -1 4 " t)," }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 51 "with boundary conditions \+ u(t, 0) = u(t,1) = 0," }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 49 "and with initial conditions u(0, x) = 0 and \+ " }{XPPEDIT 18 0 "diff(u,t);" "6#-%%diffG6$%\"uG%\"tG" }{TEXT -1 10 "( 0,x) = 0." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 69 "First, we draw the graph with Maple' s built-in numerical techniqu es." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "pde1:=diff(u(t,x),x,x )=diff(u(t,x),t,t)-16*sin(Pi*t);" }}}{PARA 0 "" 0 "" {TEXT -1 55 "Here are the boundary conditions and initial condition." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "bc1:=\{u(t,0)=0,u(t,1)=0,u(0,x)=0,D[1](u)(0 ,x)=0\};" }}}{PARA 0 "" 0 "" {TEXT -1 118 "Now, we call the numerical \+ solver for this pde. The output is a module. We recall how to make a g raph of the solution." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "u1: =pdsolve(pde1,bc1,type=numeric,time=t,range=0..1,timestep=1/100,spaces tep=1/6);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 81 "u1:-plot3d(u(t ,x),t=0..2,axes=normal,orientation=[25,65],\n style=wireframe);" }}}{PARA 0 "" 0 "" {TEXT -1 77 "Now, we construct the discrete replace ment scheme and compare with the above." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "f:=x->0:\ng:=x->0:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "F:=(t,x)->16*sin(t*Pi);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "N:=4;\nM:=4;\nk:=1/M;\ndelta:=(1/N)/(1/M);" }}}{PARA 0 "" 0 "" {TEXT -1 14 "We check that " }{XPPEDIT 18 0 "delta;" "6#%&de ltaG" }{TEXT -1 17 " is small enough." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "is(delta > 0); is(delta<=1);" }}}{PARA 0 "" 0 "" {TEXT -1 58 "Here is the matrix that is used to advance the time steps ." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 145 "A:=Matrix([[0,0,0,0,0] ,[delta^2,1-delta^2,delta^2,0,0],\n [0,delta^2,1-delta^2,delta^2,0 ],\n [0,0,delta^2,1-delta^2,delta^2],[0,0,0,0,0]]);" }}}{PARA 0 " " 0 "" {TEXT -1 40 "We need to input the initial conditions." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "U[0]:=Vector([seq(f(n/N),n=0 ..N)]);\nGINT:=Vector([seq(g(n/N),n=0..N)]);" }}}{PARA 0 "" 0 "" {TEXT -1 33 "Here is the non-homogeneous term." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "for m from 0 to 2*M do\n VF[m]:=Vector([0,seq(F (m/M,n/N),n=0..N-2),0]):\nod:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "U[1]:=A.U[0]+k^2*VF[0]/2+k*GINT;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "for m from 2 to 2*M do\n U[m]:=A.U[m-1]-U[m-2]+k^2*V F[m-1]:\nod:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 155 "for m from 1 to 2*M do J[m]:=pointplot3d([seq([(n-1)/N,m/(M),U[m][n]],n=1..N+1)] ,axes=normal,connect=true,thickness=3,color=BLACK,orientation=[25,65]) :\nod:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 87 "Jlim:=u1:-plot3d( u(t,x),t=0..2,axes=normal,orientation=[25,65],\n style=wireframe) :" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "display3d([seq(J[p],p= 1..2*M),Jlim]);" }}}{PARA 0 "" 0 "" {TEXT -1 83 "This makes a good fit to the graph Maple found with the built-in numerical methods." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}}{SECT 0 {PARA 3 "" 0 "" {TEXT 266 9 "Problem 2" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 169 "The next problem represents a string at rest and tied do wn at one end. The action is that we move the left end up and down and watch as the wave runs through the string." }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 40 "pde1:=diff(u(t,x),x,x)=diff(u(t,x),t,t);" }}}{PARA 0 "" 0 "" {TEXT -1 55 "Here are the boundary conditions and initial co ndition." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 57 "bc1:=\{u(t,0)=si n(Pi*t),u(t,1)=0,u(0,x)=0,D[1](u)(0,x)=0\};" }}}{PARA 0 "" 0 "" {TEXT -1 101 "Now, we call the numerical solver for this pde. The output is \+ a module. We recall how to make a plot." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "u1:=pdsolve(pde1,bc1,type=numeric,time=t,range=0..1,t imestep=1/100,spacestep=1/6);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 81 "u1:-plot3d(u(t,x),t=0..4,axes=normal,orientation=[25,65],\n \+ style=wireframe);" }}}{PARA 0 "" 0 "" {TEXT -1 44 "An animation reall y shows what is happening." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 103 "u1:-animate(u(t,x),t=0..5,frames=30,\n labels=[\"x\" ,\"u(x,t)\"], labelfont=[TIMES,ROMAN,14]);\n" }}}{PARA 0 "" 0 "" {TEXT -1 78 "Now, we follow methods such as we have seen to make a num erical approximation." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "f:= x->0:\ng:=x->0:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "F:=(t,x)->0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "N:=4;\nM:=4;\nk:=1/M;\ndelta:=(1/N) /(1/M);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "is(delta > 0); i s(delta<=1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 145 "A:=Matrix( [[0,0,0,0,0],[delta^2,1-delta^2,delta^2,0,0],\n [0,delta^2,1-delta ^2,delta^2,0],\n [0,0,delta^2,1-delta^2,delta^2],[0,0,0,0,0]]);" } }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "U[0]:=Vector([seq(f(n/N),n =0..N)]);\nGINT:=Vector([seq(g(n/N),n=0..N)]);" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 75 "for m from 0 to 2*M do\n VF[m]:=Vector([0,seq (F(m/M,n/N),n=0..N-2),0]):\nod:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "Vector([sin(1/4*Pi),seq(0,i=1..N)]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 80 "U[1]:=A.U[0]+k^2*VF[0]/2+k*GINT+\n Vec tor([sin(1/4*Pi),seq(0,i=1..N)]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 121 "for m from 2 to 2*M do\n U[m]:=A.U[m-1]-U[m-2]+k^2* VF[m-1]+\n Vector([sin(m/M*Pi)+U[m-2][1],seq(0,i=1..N)]):\no d:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 163 "for m from 1 to 2*M \+ do J[m]:=pointplot3d([seq([(n-1)/N,m/(M),U[m][n]],n=1..N+1)],axes=norm al,connect=true,thickness=3,color=BLACK,orientation=[25,65]):\nod:\nm: ='m':" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 87 "Jlim:=u1:-plot3d(u (t,x),t=0..2,axes=normal,orientation=[25,65],\n style=wireframe): " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "display3d([seq(J[p],p=1 ..2*M),Jlim]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {PARA 0 "" 0 "" {TEXT -1 83 "Again, this seems to be a good fit to the graph found by Maple's numerical methods." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}}{SECT 0 {PARA 3 "" 0 " " {TEXT 267 9 "Problem 3" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 50 "Consider this third partial differential equation. " }}{PARA 0 "" 0 "" {TEXT -1 3 " " }}{PARA 0 "" 0 "" {TEXT -1 8 " \+ " }{XPPEDIT 18 0 "diff(u,`$`(x,2));" "6#-%%diffG6$%\"uG-%\"$G6$% \"xG\"\"#" }{TEXT -1 3 " = " }{XPPEDIT 18 0 "diff(u,`$`(t,2));" "6#-%% diffG6$%\"uG-%\"$G6$%\"tG\"\"#" }{TEXT -1 9 " + u(t,x)" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 51 "with boundary condit ions u(t, 0) = u(t,1) = 0," }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 66 "and with initial conditions u(0, x) \+ = 1 - 2 | x - 1/2 | and " }{XPPEDIT 18 0 "diff(u,t);" "6#-%%diffG6$% \"uG%\"tG" }{TEXT -1 10 "(0,x) = 0." }}{PARA 0 "" 0 "" {TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 93 "This is a problem for which we can draw a graph of the solution with three different methods:" }}{PARA 0 "" 0 "" {TEXT -1 80 " 1. Fourier Series methods will provide an analyti c solution for this problem." }}{PARA 0 "" 0 "" {TEXT -1 116 " 2. Th e Discrete Replacement Methods will provide an approximation for the g raph of the solution for this problem." }}{PARA 0 "" 0 "" {TEXT -1 134 " 3. Maple's PDSolve routine will not only provide an approximat ion for the graph, it can also provide an animation for the solution. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 116 "Conc erning the first method, we have computed the Fourier solution. The fi rst five terms of the infinite series are" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 102 "w:=(t,x)->8/Pi^2* sum((-1)^(n+1)/(2*n-1)^2*\n cos(sqrt(1+(2*n-1)^2*Pi^2)*t)*sin((2*n-1 )*Pi*x),n=1..5);" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 102 "Concerning the second method, note that the replacement \+ equations take only a slightly different form." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 92 "eq:=(U[m](n+1)-2*U[m](n)+U[m](n-1))/H^2 =\n \+ (U[m+1](n)-2*U[m](n)+U[m-1](n))/K^2+U[m](n);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "solve(eq,U[m+1](n));" }}}{PARA 0 "" 0 "" {TEXT -1 107 "It is unlikely that programming this iteration scheme will add significantly to the ideas of this section. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 177 "However, it does seem we ll to see an animation of how the string vibrates. We use Maple's buil t in PDE solver to draw this. Already, we know how to get a graph for \+ the solution." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "pde1:=diff( u(t,x),x,x)=diff(u(t,x),t,t)+u(t,x);" }}{PARA 0 "" 0 "" {TEXT -1 55 "H ere are the boundary conditions and initial condition." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "bc1:=\{u(t,0)=0,u(t,1)=0,u(0,x)=1-2 *abs(x-1/2),D[1](u)(0,x)=0\};" }}}{PARA 0 "" 0 "" {TEXT -1 105 "Now, w e call the numerical solver for this pde. The output is a module. We i llustrate how to make a plot." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 83 "u1:=pdsolve(pde1,bc1,type=numeric,time=t,range=0..1,timestep=1/1 00,spacestep=1/20);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 59 "u1:- plot3d(u(t,x),t=0..5,axes=normal,orientation=[-35,65]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 125 "u1:-animate(u(t,x),t=0..5,frames=3 0,\n labels=[\"x\",\"u(x,t)\"], labelfont=[TIMES,ROMAN,14], \n \+ scaling=constrained);\n" }}}{PARA 0 "" 0 "" {TEXT -1 211 "It has \+ been the intent to show the basic methods for numerical methods for so lving wave equations. But also, for the user, it is good to see the in creased utility of using the numerical methods built into Maple." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 50 "EMAIL: herod@math.gatech.edu or jherod@tds.net" }}{PARA 0 " " 0 "" {TEXT -1 38 "URL: http://www.math.gatech.edu/~herod" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 256 "" 0 "" {TEXT -1 36 "Copyright \251 \+ 2003 by James V. Herod" }}{PARA 256 "" 0 "" {TEXT -1 19 "All rights \+ reserved" }}}{MARK "0 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 } {PAGENUMBERS 1 1 2 33 1 1 }