{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 "Hyperlink" -1 17 "" 0 1 0 128 128 1 2 0 1 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 14 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 257 "" 1 12 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 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" 18 264 "" 0 1 0 0 0 0 0 1 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 "MaplePi" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE " " -1 267 "MaplePi" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 268 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 269 "MaplePi" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 270 "MaplePi" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 271 "MaplePi" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 272 "MaplePi" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 273 "" 1 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 274 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 275 "MaplePi " 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 276 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 277 "" 1 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 278 "" 1 12 0 0 0 0 0 0 0 0 0 0 0 0 0 0 }{CSTYLE " " -1 279 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 282 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 283 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 284 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 285 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{PSTYLE "Normal" -1 0 1 {CSTYLE "" -1 -1 "Tim es" 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 "Norma l" -1 256 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 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Normal" -1 257 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 256 "" 0 "" {TEXT -1 0 "" }}{PARA 256 "" 0 "" {TEXT 256 88 "Sect ion 8.1: Numerical Solutions for Ordinary Differential Equations: Diff erence Methods" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 0 {PARA 3 "" 0 "" {TEXT 257 34 "Maple Packages used in Section 8.1" }}{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(LinearAlgebra):" }}}}{PARA 0 "" 0 "" {TEXT -1 1 " " }} {PARA 0 "" 0 "" {TEXT -1 474 "There are many methods for getting numer ical solutions for ordinary differential equations. Through the years, with improving analysis, software and hardware, the methods have impr oved in accuracy. Euler's Method, Improved Euler's Method, and the Run ge-Kutta Method are commonly taught in undergraduate classes as method s for obtaining numerical solutions for ordinary differential equation s. Maple uses a Runge-Kutta-Fehlberg 4/5th order integration method, \+ among others. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 187 "In this Section 8.1, we will follow closely the methods \+ presented by Douglas Meade in Unit 18 of his notes for the Ordinary Di fferential Equations Powertool at the Maplesoft web site. See" }} {PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {URLLINK 17 "http://ww w.mapleapps.com/powertools/des/des.shtml" 4 "" "" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 191 "In that material, Meade \+ has a complete section, Part IV, on numerical methods for initial valu e value problems. Part V in his materials contains numerical methods f or boundary value problems." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 232 "Of course, in partial differential equations, \+ derivatives are made with respect to more than one variable, so the ab ility to simply \"step across the interval\" is missing. The methods s hown in this 8th Chapter uses the techniques of " }{TEXT 258 25 "finit e difference methods" }{TEXT -1 30 ". The methods could be called " } {TEXT 259 28 "discrete replacement methods" }{TEXT -1 63 " for, follow ing the development in Unit 18 of Meade's work, we " }{TEXT 260 7 "rep lace" }{TEXT -1 20 " the derivatives by " }{TEXT 276 8 "discrete" } {TEXT -1 154 " approximations. With these methods, the job of getting \+ an approximation for a solution of a differential equation is changed \+ to a task in linear algebra." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 151 "The purpose of reviewing these ordinary differ ential equations methods is that similar techniques will be expanded f or partial differential equations. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" } }{PARA 0 "" 0 "" {TEXT -1 345 "It should be no surprise that there are choices on what discrete replacements will be made for the derivative s. For example, there are forward difference replacements and backward difference replacements. We follow Meade as we make choices for the r eplacements. We use and apply the finite difference method to a two-po int boundary value problem" }}{PARA 257 "" 0 "" {TEXT -1 2 " " } {XPPEDIT 18 0 "Diff(y,x$2) + p(x)*Diff(y,x) + q(x)*y = f(x)" "6#/,(-%% DiffG6$%\"yG-%\"$G6$%\"xG\"\"#\"\"\"*&-%\"pG6#F,F.-F&6$F(F,F.F.*&-%\"q G6#F,F.F(F.F.-%\"fG6#F," }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 127 "There will be a boundary condition \+ at each endpoint of an interval over which we solve the equation. Mead e calls the endpoints " }{TEXT 262 1 "a" }{TEXT -1 5 " and " }{TEXT 261 1 "b" }{TEXT -1 35 ". The mesh points will be denoted [" } {XPPEDIT 18 0 "x[0];" "6#&%\"xG6#\"\"!" }{TEXT -1 2 ", " }{XPPEDIT 18 0 "x[1],x[2];" "6$&%\"xG6#\"\"\"&F$6#\"\"#" }{TEXT -1 8 ", ... , " } {XPPEDIT 18 0 "x[N];" "6#&%\"xG6#%\"NG" }{TEXT -1 2 "]." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 35 "We explain now how w e will replace " }{XPPEDIT 18 0 "diff(y,x);" "6#-%%diffG6$%\"yG%\"xG" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "diff(y,`$`(x,2));" "6#-%%diffG6$% \"yG-%\"$G6$%\"xG\"\"#" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT 263 15 "Replacement for" }{TEXT -1 2 " " }{XPPEDIT 264 0 "diff(y,x);" "6# -%%diffG6$%\"yG%\"xG" }{TEXT -1 1 ":" }}{PARA 0 "" 0 "" {TEXT -1 23 "R ecall Taylor's Series:" }}{PARA 0 "" 0 "" {TEXT -1 22 " y( x + h) = \+ y(x) + " }{TEXT 265 1 "h" }{TEXT -1 11 " y '(x) + " }{TEXT 266 1 "0" }{TEXT -1 2 "( " }{XPPEDIT 18 0 "h^2;" "6#*$%\"hG\"\"#" }{TEXT -1 1 ") " }}{PARA 0 "" 0 "" {TEXT -1 3 "and" }}{PARA 0 "" 0 "" {TEXT -1 22 " \+ y( x - h) = y(x) - " }{TEXT 274 1 "h" }{TEXT -1 11 " y '(x) + " } {TEXT 275 1 "0" }{TEXT -1 2 "( " }{XPPEDIT 18 0 "h^2;" "6#*$%\"hG\"\"# " }{TEXT -1 1 ")" }}{PARA 0 "" 0 "" {TEXT -1 5 "Thus," }}{PARA 0 "" 0 "" {TEXT -1 40 " y '(x) is replaced by the difference " }{XPPEDIT 18 0 "1/2;" "6#*&\"\"\"F$\"\"#!\"\"" }{TEXT -1 1 " " }{XPPEDIT 18 0 "( y(x+h)-y(x-h))/h;" "6#*&,&-%\"yG6#,&%\"xG\"\"\"%\"hGF*F*-F&6#,&F)F*F+! \"\"F/F*F+F/" }{TEXT -1 31 " and this is good with order " }{TEXT 267 1 "0" }{TEXT -1 2 "( " }{XPPEDIT 18 0 "h^2;" "6#*$%\"hG\"\"#" } {TEXT -1 2 ")." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT 268 15 "Replacement for" }{TEXT -1 2 " " }{XPPEDIT 257 0 "diff( y,`$`(x,2));" "6#-%%diffG6$%\"yG-%\"$G6$%\"xG\"\"#" }{TEXT -1 1 ":" }} {PARA 0 "" 0 "" {TEXT -1 26 "Use Taylor's Series again:" }}{PARA 0 "" 0 "" {TEXT -1 32 " y(x + h) = y(x) + h y '(x) + " }{XPPEDIT 18 0 "h^ 2/2;" "6#*&%\"hG\"\"#F%!\"\"" }{TEXT -1 11 " y ''(x) + " }{XPPEDIT 18 0 "h^3/3!;" "6#*&%\"hG\"\"$-%*factorialG6#F%!\"\"" }{TEXT -1 12 " y '' '(x) + " }{TEXT 269 1 "0" }{TEXT -1 2 "( " }{XPPEDIT 18 0 "h^4;" "6#*$ %\"hG\"\"%" }{TEXT -1 35 ").\n y(x - h) = y(x) - h y '(x) + " } {XPPEDIT 18 0 "h^2/2;" "6#*&%\"hG\"\"#F%!\"\"" }{TEXT -1 12 " y ''(x) \+ - " }{XPPEDIT 18 0 "h^3/3!;" "6#*&%\"hG\"\"$-%*factorialG6#F%!\"\"" } {TEXT -1 12 " y '''(x) + " }{TEXT 270 1 "0" }{TEXT -1 2 "( " } {XPPEDIT 18 0 "h^4;" "6#*$%\"hG\"\"%" }{TEXT -1 2 ")." }}{PARA 0 "" 0 "" {TEXT -1 37 "Add these two equations and we we get" }}{PARA 0 "" 0 "" {TEXT -1 41 " y( x + h) - y(x) + y(x - h) - y(x) = " }{XPPEDIT 18 0 "h^2;" "6#*$%\"hG\"\"#" }{TEXT -1 12 " y '' (x) + " }{TEXT 271 1 "0" }{TEXT -1 2 "( " }{XPPEDIT 18 0 "h^4;" "6#*$%\"hG\"\"%" }{TEXT -1 2 ")." }}{PARA 0 "" 0 "" {TEXT -1 5 "Thus," }}{PARA 0 "" 0 "" {TEXT -1 43 " y ''(x) is replaced by the difference " }{XPPEDIT 18 0 "(y (x+h)-2*y(x)+y(x-h))/(h^2);" "6#*&,(-%\"yG6#,&%\"xG\"\"\"%\"hGF*F**&\" \"#F*-F&6#F)F*!\"\"-F&6#,&F)F*F+F0F*F**$F+F-F0" }{TEXT -1 30 " and th is is good with order " }{TEXT 272 1 "0" }{TEXT -1 2 "( " }{XPPEDIT 18 0 "h^2;" "6#*$%\"hG\"\"#" }{TEXT -1 2 ")." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 247 "We now work three problems. Th e first two were worked in Meade's Unit 18. We follow his syntax almos t in the manner of \"cut and paste.\" In this Section 8.1, before the problems are finished we will have restated each of them as a matrix \+ problem. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 0 {PARA 3 "" 0 "" {TEXT 273 9 "Example 1" }}{PARA 0 "" 0 "" {TEXT -1 35 "Consider the bo undary value problem" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 77 "ode1 :=diff(y(x),x$2)+2*diff(y(x),x)-10*y(x)=7*exp(-x)+4;\nbc1:=y(0)=2,y(1) =-5;" }}}{PARA 0 "" 0 "" {TEXT -1 48 "with boundary conditions specifi ed at the points" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "a:=0;\nb :=1;" }}}{PARA 0 "" 0 "" {TEXT -1 275 "To find the finite difference s olution associated with this boundary value problem, Meade chooses 11 \+ mesh points from a to b. He replaces the derivatives as we indicated a bove and shows graphically what a good approximation this makes to the exact solution. See his Unit 18. " }}{PARA 0 "" 0 "" {TEXT -1 229 "Ra ther than repeat what he has done, we choose only 6 mesh points and em phasize the structure of the resulting equations to be solved. We will find there is a tridiagonal matrix associated with the analysis. We p roceed. Choose N." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "N:=5;" } }}{PARA 0 "" 0 "" {TEXT -1 52 "This defines h and enables us to comput e the X mesh." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "h:=(b-a)/N; \nX:=k->a+k*h;" }}}{PARA 0 "" 0 "" {TEXT -1 56 "Next, we define the fi nite difference replacement terms." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 63 "Yp:=k->(y[k+1]-y[k-1])/2/h;\nYpp:=k->(y[k+1]-2*y[k]+y [k-1])/h^2;" }}}{PARA 0 "" 0 "" {TEXT -1 46 "The boundary conditions p rovide two equations." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "eq[ 0]:=y[0]=rhs(bc1[1]);\neq[N]:=y[N]= rhs(bc1[2]);" }}}{PARA 0 "" 0 "" {TEXT -1 118 "Here are the remaining N-1 equations for the values at t he interior mesh points. We make these equations by replacing " }} {PARA 0 "" 0 "" {TEXT -1 48 " y '(x) and y ''(x)." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 117 "for k from 1 to \+ N-1 do\n eq[k]:=eval(ode1,\{x=X(k),y(x)=y[k],\n diff(y(x),x)= Yp(k),diff(y(x),x$2)=Ypp(k)\});\nod;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 167 "This system of N+1 \+ equations can be realized in a tridiagonal matrix equation. The matrix will need N+1 rows. The first one comes from the boundary condition a t x = a." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "s[0]:=[1,seq(0,p =1..N)];" }}}{PARA 0 "" 0 "" {TEXT -1 53 "The next N-1 rows come from \+ the interior mesh points." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 149 "for n from 1 to N-1 do\n s[n]:=[seq(0,p=1..n-1),\ncoeff(lhs(eq[ n]),y[n-1]),\ncoeff(lhs(eq[n]),y[n]),\ncoeff(lhs(eq[n]),y[n+1]),\nseq( 0,p=1..N-1-n)];\nod;" }}}{PARA 0 "" 0 "" {TEXT -1 79 "Finally, there i s the last row that comes from the boundary condition at x = b." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "s[N]:=[seq(0,p=1..N),1];" }} }{PARA 0 "" 0 "" {TEXT -1 62 "Now, we write this as a matrix. We list \+ out the rows together." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "ro se:=[s[0],seq(s[n],n=1..N-1),s[N]];" }}}{PARA 0 "" 0 "" {TEXT -1 77 "W e construct the matrix with these rows. Note that A is a tridiagonal m atrix." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "A:=Matrix(rose);" }}}{PARA 0 "" 0 "" {TEXT -1 28 "Here are the unknown values " } {XPPEDIT 18 0 "y[1],y[2];" "6$&%\"yG6#\"\"\"&F$6#\"\"#" }{TEXT -1 7 ", ... ," }{XPPEDIT 18 0 "y[N];" "6#&%\"yG6#%\"NG" }{TEXT -1 1 "." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "Y:=Vector([seq(y[p],p=0..N)] );" }}}{PARA 0 "" 0 "" {TEXT -1 72 "Here is the right hand side of the equations written as a column vector." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "V:=Vector([rhs(bc1[1]),seq(rhs(eq[p]),p=1..N-1),rhs(b c1[2])]);" }}}{PARA 0 "" 0 "" {TEXT -1 48 "Finally, here is the equati on to be solve for y." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "A.Y= V;" }}}{PARA 0 "" 0 "" {TEXT -1 73 "There are many ways to solve matri x equations. We leave this up to Maple." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Y:=LinearSolve(A,V);" }}}{PARA 0 "" 0 "" {TEXT -1 80 "In an attempt to make these numbers understandable by humans, we give a listing." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "result:=eval( [seq([X(k-1),evalf(Y[k],5)],k=1..N+1)]);" }}}{PARA 0 "" 0 "" {TEXT -1 230 "It will be interesting to compare how good this approximation is \+ with the exact solution. The approximation may be improved by rerunnin g this example and increasing the size of N.\nWe draw graphs. First, c ompute the exact solution." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "sol1:=combine(dsolve(\{ode1,bc1\},y(x)));" }}}{PARA 0 "" 0 "" {TEXT -1 64 "Next, we superimpose the exact solution with this approxi mation." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 111 "plot([rhs(sol1), result],x=a..b,\n color=[BLACK,GREEN],thickness=[1,3],\n legend=[`ex act`,`finite difference`]);" }}}{PARA 0 "" 0 "" {TEXT -1 108 "Now a ba d fit, yes? I say again, you might go back and increase the size of N \+ to improve the approximation. " }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {SECT 0 {PARA 3 "" 0 "" {TEXT 277 9 "Example 2" }}{PARA 0 "" 0 "" {TEXT -1 111 "Here is the second example from Meade's Unit 18. We make modifications to see it as a tridiagonal matrix again." }}{PARA 0 "" 0 "" {TEXT -1 79 "The boundary value problem in this example does not \+ have constant coefficients." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 54 "ode2:=x^2*diff(y(x),x$2)+x*diff(y(x),x)+4*y(x)=-2*x+7;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "bc2:=y(1)=7,y(4)=-1;" }}}{PARA 0 "" 0 "" {TEXT -1 36 "The boundary conditions are given at" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "a:=1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "b:= 4;" }}}{PARA 0 "" 0 "" {TEXT -1 89 "Here is where we choose the number of mesh points. The number of mesh points will be N+1." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "N:=5;" }}}{PARA 0 "" 0 "" {TEXT -1 116 "The interval over which we approximate this equation is [a, b] an d for this interval, we have defined the step size." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "h:=(b-a)/N;" }}}{PARA 0 "" 0 "" {TEXT -1 25 "Here are the mesh points." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "X:=k->a+k*h;" }}}{PARA 0 "" 0 "" {TEXT -1 77 "The replacement op erators are as they were explained in the discussion above." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "Yp:=k->(y[k+1]-y[k-1])/(2*h);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Ypp:=k->(y[k+1]-2*y[k]+y[k-1 ])/h^2;" }}}{PARA 0 "" 0 "" {TEXT -1 46 "The boundary conditions provi de two equations." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "eq[0]:= y[0]=rhs(bc2[1]);\neq[N]:=y[N]= rhs(bc2[2]);" }}}{PARA 0 "" 0 "" {TEXT -1 118 "Here are the remaining N-1 equations for the values at t he interior mesh points. We make these equations by replacing " }} {PARA 0 "" 0 "" {TEXT -1 48 " y '(x) and y ''(x)." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 120 "for k from 1 to \+ N-1 do\n eq[k]:=eval(ode2,\{x=X(k),y(x)=y[k],\n diff(y(x),x )=Yp(k),diff(y(x),x$2)=Ypp(k)\} );\nod;" }}}{PARA 0 "" 0 "" {TEXT -1 160 "We realize this system of N+1 equations as tridiagonal matrix equ ation. The matrix will need N+1 rows. The first one comes from the bou ndary condition at x = a." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "s[0]:=[1,seq(0,p=1..N)];" }}}{PARA 0 "" 0 "" {TEXT -1 53 "The next N- 1 rows come from the interior mesh points." }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 149 "for n from 1 to N-1 do\n s[n]:=[seq(0,p=1..n-1), \ncoeff(lhs(eq[n]),y[n-1]),\ncoeff(lhs(eq[n]),y[n]),\ncoeff(lhs(eq[n]) ,y[n+1]),\nseq(0,p=1..N-1-n)];\nod;" }}}{PARA 0 "" 0 "" {TEXT -1 79 "F inally, there is the last row that comes from the boundary condition a t x = b." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "s[N]:=[seq(0,p=1 ..N),1];" }}}{PARA 0 "" 0 "" {TEXT -1 62 "Now, we write this as a matr ix. We list out the rows together." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "rose:=[s[0],seq(s[n],n=1..N-1),s[N]];" }}}{PARA 0 "" 0 "" {TEXT -1 40 "We construct the matrix with these rows." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "A:=Matrix(rose);" }}}{PARA 0 "" 0 " " {TEXT -1 28 "Here are the unknown values " }{XPPEDIT 18 0 "y[1],y[2] ;" "6$&%\"yG6#\"\"\"&F$6#\"\"#" }{TEXT -1 7 ", ... ," }{XPPEDIT 18 0 " y[N];" "6#&%\"yG6#%\"NG" }{TEXT -1 1 "." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "Y:=Vector([seq(y[p],p=0..N)]);" }}}{PARA 0 "" 0 "" {TEXT -1 72 "Here is the right hand side of the equations written as a column vector." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "V:=Vector ([rhs(bc2[1]),seq(rhs(eq[p]),p=1..N-1),rhs(bc2[2])]);" }}}{PARA 0 "" 0 "" {TEXT -1 48 "Finally, here is the equation to be solve for y." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "A.Y=V;" }}}{PARA 0 "" 0 "" {TEXT -1 36 "Maple will solve this system for us." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Y:=LinearSolve(A,V);" }{TEXT -1 0 "" }}} {PARA 0 "" 0 "" {TEXT -1 22 "We reprint for humans." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "result:=eval([seq([X(k-1),evalf(Y[k],5)],k= 1..N+1)]);" }}}{PARA 0 "" 0 "" {TEXT -1 60 "\nFinally, we draw graphs. First, compute the exact solution." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "sol2:=combine(dsolve(\{ode2,bc2\},y(x)));" }}}{PARA 0 "" 0 "" {TEXT -1 64 "Next, we superimpose the exact solution with th is approximation." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 114 "plot([ rhs(sol2),result],x=a..b,\n color=[BLACK,GREEN], thickness=[1,3],\n \+ legend=[`exact`,`finite difference`]);" }}}{PARA 0 "" 0 "" {TEXT -1 102 "The fit with N = 5 is not as good as we would have liked. Going b ack and increasing N is irresistible." }}{PARA 0 "" 0 "" {TEXT -1 0 " " }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 0 {PARA 3 "" 0 "" {TEXT 278 9 "Example 3" }}{PARA 0 "" 0 "" {TEXT -1 220 "This example differs from the previous two in that the boundary condition at the right han d side involves the derivative of y, instead of the value of y. Also, \+ we make the right hand side of the equation not continuous. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 71 "Here is the exa mple: we seek a graph for a function y(x) that satisfies" }}{PARA 0 " " 0 "" {TEXT -1 11 " " }{XPPEDIT 18 0 "diff(y(x),`$`(x,2)); " "6#-%%diffG6$-%\"yG6#%\"xG-%\"$G6$F)\"\"#" }{TEXT -1 16 " - 3 y(x) = f(x)" }}{PARA 0 "" 0 "" {TEXT -1 4 "with" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 81 "boundary conditions y(-1) = 1 and \+ y '(1) = -1. The function f is indicated below." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 75 "To the extent possible, w e follow the pattern of the previous two examples." }}{PARA 0 "" 0 "" {TEXT -1 68 "We use a second order differential equation with forcing \+ function f." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "ode3:=diff(y( x),x$2)-9*y(x)=f(x);\nbc3:=y(-1)=1,D(y)(1)=-1;" }}}{EXCHG {PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 22 "f:=x->(1-signum(x))/2;" }}}{PARA 0 "" 0 "" {TEXT -1 36 "The boundary conditions are given at" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "a:=-1;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "b: =1;" }}}{PARA 0 "" 0 "" {TEXT -1 89 "Here is where we choose the numbe r of mesh points. The number of mesh points will be N+1." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "N:=9;" }}}{PARA 0 "" 0 "" {TEXT -1 116 "The interval over which we approximate this equation is [a, b] an d for this interval, we have defined the step size." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "h:=(b-a)/N;" }}}{PARA 0 "" 0 "" {TEXT -1 25 "Here are the mesh points." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "X:=k->a+k*h;" }}}{PARA 0 "" 0 "" {TEXT -1 77 "The replacement op erators are as they were explained in the discussion above." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "Yp:=k->(y[k+1]-y[k-1])/(2*h);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "Ypp:=k->(y[k+1]-2*y[k]+y[k-1 ])/h^2;" }}}{PARA 0 "" 0 "" {TEXT -1 37 "Of course, the boundary condi tion at " }{TEXT 279 1 "a" }{TEXT -1 103 " provides an equation of the same character as in the previous two examples. The boundary conditio n at " }{TEXT 280 1 "b" }{TEXT -1 23 " needs some discussion." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "eq[0]:=y[0]=rhs(bc3[1]);" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 95 "Here, we explain how to get equ ation N from the boundary condition at the right hand endpoint, " } {TEXT 281 1 "b" }{TEXT -1 47 ". Our replacement scheme for the derivat ive at " }{TEXT 282 1 "b" }{TEXT -1 7 " = 1 = " }{XPPEDIT 18 0 "X[N]; " "6#&%\"XG6#%\"NG" }{TEXT -1 11 " should be" }}{PARA 0 "" 0 "" {TEXT -1 17 " " }{XPPEDIT 18 0 "(Y[N+1]-Y[N-1])/(2*h); " "6#*&,&&%\"YG6#,&%\"NG\"\"\"F*F*F*&F&6#,&F)F*F*!\"\"F.F**&\"\"#F*%\" hGF*F." }{TEXT -1 8 " = -1." }}{PARA 0 "" 0 "" {TEXT -1 25 "However, the sequence of " }{TEXT 283 2 "Y " }{TEXT -1 18 "'s runs only from \+ " }{XPPEDIT 18 0 "Y[0];" "6#&%\"YG6#\"\"!" }{TEXT -1 4 " to " } {XPPEDIT 18 0 "Y[N];" "6#&%\"YG6#%\"NG" }{TEXT -1 15 " . Thus we get \+ " }{XPPEDIT 18 0 "eq[N];" "6#&%#eqG6#%\"NG" }{TEXT -1 61 " in the foll owing manner. From the previous line, we get that" }}{PARA 0 "" 0 "" {TEXT -1 10 " " }{XPPEDIT 18 0 "Y[N+1];" "6#&%\"YG6#,&%\"NG\" \"\"F(F(" }{TEXT -1 3 " = " }{XPPEDIT 18 0 "Y[n-1];" "6#&%\"YG6#,&%\"n G\"\"\"F(!\"\"" }{TEXT -1 5 " - 2 " }{TEXT 284 1 "h" }{TEXT -1 1 "." } }{PARA 0 "" 0 "" {TEXT -1 50 "If we followed the discrete replacement \+ scheme at " }{XPPEDIT 18 0 "X[N];" "6#&%\"XG6#%\"NG" }{TEXT -1 15 ", w e would have" }}{PARA 0 "" 0 "" {TEXT -1 7 " " }{XPPEDIT 18 0 "( Y[N+1]-2*Y[N]+Y[N-1])/(h^2);" "6#*&,(&%\"YG6#,&%\"NG\"\"\"F*F*F**&\"\" #F*&F&6#F)F*!\"\"&F&6#,&F)F*F*F/F*F**$%\"hGF,F/" }{TEXT -1 6 " - 9 " }{XPPEDIT 18 0 "Y[N];" "6#&%\"YG6#%\"NG" }{TEXT -1 7 " = f( " } {XPPEDIT 18 0 "X[N];" "6#&%\"XG6#%\"NG" }{TEXT -1 3 " )." }}{PARA 0 " " 0 "" {TEXT -1 27 "Again, this has a value of " }{TEXT 285 1 "Y" } {TEXT -1 9 ", namely " }{XPPEDIT 18 0 "Y[N+1];" "6#&%\"YG6#,&%\"NG\"\" \"F(F(" }{TEXT -1 53 ", that we do not use. It is replaced by the valu e of " }{XPPEDIT 18 0 "Y[N+1];" "6#&%\"YG6#,&%\"NG\"\"\"F(F(" }{TEXT -1 49 " that we got from the boundary condition to yield" }}{PARA 0 " " 0 "" {TEXT -1 6 " " }{XPPEDIT 18 0 "(Y[N-1]-2*h-2*Y[N]+Y[N-1])/ (h^2);" "6#*&,*&%\"YG6#,&%\"NG\"\"\"F*!\"\"F**&\"\"#F*%\"hGF*F+*&F-F*& F&6#F)F*F+&F&6#,&F)F*F*F+F*F**$F.F-F+" }{TEXT -1 6 " - 9 " }{XPPEDIT 18 0 "Y[N];" "6#&%\"YG6#%\"NG" }{TEXT -1 7 " = f( " }{XPPEDIT 18 0 "X [N];" "6#&%\"XG6#%\"NG" }{TEXT -1 17 " ).\nThis defines " }{XPPEDIT 18 0 "eq[N];" "6#&%#eqG6#%\"NG" }{TEXT -1 1 "." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 71 "eq[N]:=simplify((y[N-1]-2*y[N]+y[N-1])/h^2-9*y[N ]=\n f(X(N))+2/h);" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 " " 0 "" {TEXT -1 118 "Here are the remaining N-1 equations for the valu es at the interior mesh points. We make these equations by replacing \+ " }}{PARA 0 "" 0 "" {TEXT -1 48 " y '(x) a nd y ''(x)." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 120 "for k from 1 to N-1 do\n eq[k]:=eval(ode3,\{x=X(k),y(x)=y[k],\n diff(y( x),x)=Yp(k),diff(y(x),x$2)=Ypp(k)\} );\nod;" }}}{PARA 0 "" 0 "" {TEXT -1 160 "We realize this system of N+1 equations as tridiagonal matrix \+ equation. The matrix will need N+1 rows. The first one comes from the \+ boundary condition at x = a." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "s[0]:=[1,seq(0,p=1..N)];" }}}{PARA 0 "" 0 "" {TEXT -1 53 "The next N-1 rows come from the interior mesh points." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 149 "for n from 1 to N-1 do\n s[n]:=[seq(0,p=1..n- 1),\ncoeff(lhs(eq[n]),y[n-1]),\ncoeff(lhs(eq[n]),y[n]),\ncoeff(lhs(eq[ n]),y[n+1]),\nseq(0,p=1..N-1-n)];\nod;" }}}{PARA 0 "" 0 "" {TEXT -1 79 "Finally, there is the last row that comes from the boundary condit ion at x = b." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 106 "s[N]:=[seq (0,p=1..N-2),coeff(lhs(eq[N]),y[N-2]),\n coeff(lhs(eq[N]),y[N-1 ]),coeff(lhs(eq[N]),y[N])];" }}}{PARA 0 "" 0 "" {TEXT -1 62 "Now, we w rite this as a matrix. We list out the rows together." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "rose:=[s[0],seq(s[n],n=1..N-1),s[N]];" }} }{PARA 0 "" 0 "" {TEXT -1 40 "We construct the matrix with these rows. " }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "A:=Matrix(rose);" }}} {PARA 0 "" 0 "" {TEXT -1 28 "Here are the unknown values " }{XPPEDIT 18 0 "y[1],y[2];" "6$&%\"yG6#\"\"\"&F$6#\"\"#" }{TEXT -1 7 ", ... ," } {XPPEDIT 18 0 "y[N];" "6#&%\"yG6#%\"NG" }{TEXT -1 1 "." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "Y:=Vector([seq(y[p],p=0..N)]);" }}} {PARA 0 "" 0 "" {TEXT -1 72 "Here is the right hand side of the equati ons written as a column vector." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 61 "V:=Vector([rhs(bc3[1]),seq(rhs(eq[p]),p=1..N-1),rhs(eq[N])]); " }}}{PARA 0 "" 0 "" {TEXT -1 48 "Finally, here is the equation to be \+ solve for y." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 6 "A.Y=V;" }}} {PARA 0 "" 0 "" {TEXT -1 36 "Maple will solve this system for us." }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "Y:=LinearSolve(A,V);" } {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 22 "We reprint for humans." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 53 "result:=eval([seq([X(k-1), evalf(Y[k],5)],k=1..N+1)]);" }}}{PARA 0 "" 0 "" {TEXT -1 60 "\nFinally , we draw graphs. First, compute the exact solution." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "sol3:=combine(dsolve(\{ode3,bc3\},y(x))); " }}}{PARA 0 "" 0 "" {TEXT -1 64 "Next, we superimpose the exact solut ion with this approximation." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 115 "plot([rhs(sol3),result],x=-1..1,\n color=[BLACK,GREEN], thickne ss=[1,3],\n legend=[`exact`,`finite difference`]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 81 "If you d on't think this fit is pretty good, go back and make N larger. Try N = 9." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 226 "This review of the finite difference m ethods, or discrete replacement methods, for boundary value problems i n differential equations was made to suggest the techniques for partia l differential equations. We examine these next." }}{PARA 0 "" 0 "" {TEXT -1 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 257 "" 0 "" {TEXT -1 36 "Copyright \251 2003 by James V. Herod" }} {PARA 257 "" 0 "" {TEXT -1 19 "All rights reserved" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{MARK "0 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 } {PAGENUMBERS 1 1 2 33 1 1 }