{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 "Normal" -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 {PARA 256 "" 0 "" {TEXT 256 88 "Section 8.1: Numerical Solutio
ns for Ordinary Differential Equations: Difference Methods" }}{PARA 0
"" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 "" {TEXT 257 34 "Maple Pa
ckages 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 475 "There \+
are many methods for getting numerical solutions for ordinary differen
tial equations. Through the years, with improving analysis, software, \+
and hardware, the methods have improved in accuracy. Euler's Method, I
mproved Euler's Method, and the Runge-Kutta Method are commonly taught
in undergraduate classes as methods for obtaining numerical solutions
for ordinary differential equations. Maple uses a Runge-Kutta-Fehlbe
rg 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 Differential Equations Powertool at the \+
Maplesoft web site. See" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 ""
0 "" {URLLINK 17 "http://www.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 numer
ical methods for initial value value problems. Part V in his materials
contains numerical methods for boundary value problems." }}{PARA 0 "
" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 232 "Of course, in par
tial differential equations, derivatives are made with respect to more
than one variable, so the ability to simply \"step across the interva
l\" is missing. The methods shown in this 8th Chapter uses the techniq
ues of " }{TEXT 258 25 "finite difference methods" }{TEXT -1 30 ". The
methods could be called " }{TEXT 259 28 "discrete replacement methods
" }{TEXT -1 63 " for, following the development in Unit 18 of Meade's \+
work, we " }{TEXT 260 7 "replace" }{TEXT -1 20 " the derivatives by "
}{TEXT 276 8 "discrete" }{TEXT -1 154 " approximations. With these met
hods, the job of getting an approximation for a solution of a differen
tial equation is changed to a task in linear algebra." }}{PARA 0 "" 0
"" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 151 "The purpose of review
ing these ordinary differential equations methods is that similar tech
niques will be expanded for partial differential equations. " }}{PARA
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 345 "It should be n
o surprise that there are choices on what discrete replacements will b
e made for the derivatives. For example, there are forward difference \+
replacements and backward difference replacements. We follow Meade as \+
we make choices for the replacements. We use and apply the finite diff
erence method to a two-point 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.*&-%\"qG6#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 sol
ve the equation. Meade 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 explai
n now how we 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 "Recall 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 or
der " }{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 this 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 w
ork three problems. The first two were worked in Meade's Unit 18. We f
ollow his syntax almost in the manner of \"cut and paste.\" In this S
ection 8.1, before the problems are finished we will have restated eac
h of them as a matrix problem. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{SECT 1 {PARA 3 "" 0 "" {TEXT 273 9 "Example 1" }}{PARA 0 "" 0 ""
{TEXT -1 35 "Consider the boundary 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 specified 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 solution associated with this boundary valu
e problem, Meade chooses 11 mesh points from a to b. He replaces the d
erivatives as we indicated above and shows graphically what a good app
roximation this makes to the exact solution. See his Unit 18. " }}
{PARA 0 "" 0 "" {TEXT -1 229 "Rather than repeat what he has done, we \+
choose only 6 mesh points and emphasize the structure of the resulting
equations to be solved. We will find there is a tridiagonal matrix as
sociated with the analysis. We proceed. Choose N." }}{EXCHG {PARA 0 ">
" 0 "" {MPLTEXT 1 0 5 "N:=5;" }}}{PARA 0 "" 0 "" {TEXT -1 52 "This de
fines h and enables us to compute 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 finite 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 provide two equations." }}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 50 "eq[0]:=y[0]=rhs(bc1[1]);\neq[N]:=y[N]= rhs(b
c1[2]);" }}}{PARA 0 "" 0 "" {TEXT -1 118 "Here are the remaining N-1 e
quations for the values at the interior mesh points. We make these equ
ations 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 fr
om 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 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 "rose:=[s[0],seq(s[n],n=1..N-1),s[N]
];" }}}{PARA 0 "" 0 "" {TEXT -1 77 "We construct the matrix with these
rows. Note that A is a tridiagonal matrix." }}{EXCHG {PARA 0 "> " 0 "
" {MPLTEXT 1 0 16 "A:=Matrix(rose);" }}}{PARA 0 "" 0 "" {TEXT -1 28 "H
ere are the unknown values " }{XPPEDIT 18 0 "y[1],y[2];" "6$&%\"yG6#\"
\"\"&F$6#\"\"#" }{TEXT -1 7 ", ... ," }{XPPEDIT 18 0 "y[N];" "6#&%\"yG
6#%\"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(r
hs(eq[p]),p=1..N-1),rhs(bc1[2])]);" }}}{PARA 0 "" 0 "" {TEXT -1 48 "Fi
nally, 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 73 "There are
many ways to solve matrix 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 under
standable 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 h
ow good this approximation is with the exact solution. The approximati
on may be improved by rerunning this example and increasing the size o
f N.\nWe draw graphs. First, compute 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 approximation." }}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 111 "plot([rhs(sol1),result],x=a..b,\n color=[BLACK,GREE
N],thickness=[1,3],\n legend=[`exact`,`finite difference`]);" }}}
{PARA 0 "" 0 "" {TEXT -1 108 "Now a bad fit, yes? I say again, you mig
ht go back and increase the size of N to improve the approximation. "
}}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{SECT 1 {PARA 3 "" 0 "" {TEXT 277
9 "Example 2" }}{PARA 0 "" 0 "" {TEXT -1 111 "Here is the second examp
le from Meade's Unit 18. We make modifications to see it as a tridiago
nal matrix again." }}{PARA 0 "" 0 "" {TEXT -1 79 "The boundary value p
roblem 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 point
s 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 t
his equation is [a, b] and 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 operators are as they were explained in t
he 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 provide two equations." }}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 50 "eq[0]:=y[0]=rhs(bc2[1]);\neq[N]:=y[N]= rhs(b
c2[2]);" }}}{PARA 0 "" 0 "" {TEXT -1 118 "Here are the remaining N-1 e
quations for the values at the interior mesh points. We make these equ
ations 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 equation. The matrix will need N+1 rows. The fir
st 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 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 "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:=Matr
ix(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 th
e 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(b
c2[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 36 "Maple will solve this system for u
s." }}{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 "\nFinal
ly, 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 this approximation." }}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 114 "plot([rhs(sol2),result],x=a..b,\n color=[BLACK,GRE
EN], 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 back and increasing N is irresistible." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{SECT 1 {PARA 3 "" 0 "" {TEXT 278 9 "Example 3" }}{PARA 0 "" 0 ""
{TEXT -1 220 "This example differs from the previous two in that the b
oundary condition at the right hand side involves the derivative of y,
instead of the value of y. Also, we make the right hand side of the e
quation not continuous. " }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "
" 0 "" {TEXT -1 71 "Here is the example: we seek a graph for a functio
n 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 i
ndicated below." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 ""
{TEXT -1 75 "To the extent possible, we follow the pattern of the prev
ious two examples." }}{PARA 0 "" 0 "" {TEXT -1 68 "We use a second ord
er 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 number of mesh points. The number of m
esh 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 appro
ximate this equation is [a, b] and 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 operators 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 condition at " }{TEXT 279 1 "a" }{TEXT
-1 103 " provides an equation of the same character as in the previous
two examples. The boundary condition 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 equation N from the boundary condition at
the right hand endpoint, " }{TEXT 281 1 "b" }{TEXT -1 47 ". Our repla
cement scheme for the derivative 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 following 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#,&%\"nG\"\"\"F(!\"\"" }{TEXT -1 5 " - 2 " }{TEXT 284 1 "h
" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 50 "If we followed the di
screte replacement scheme at " }{XPPEDIT 18 0 "X[N];" "6#&%\"XG6#%\"NG
" }{TEXT -1 15 ", we would have" }}{PARA 0 "" 0 "" {TEXT -1 7 " \+
" }{XPPEDIT 18 0 "(Y[N+1]-2*Y[N]+Y[N-1])/(h^2);" "6#*&,(&%\"YG6#,&%\"N
G\"\"\"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 rep
laced by the value of " }{XPPEDIT 18 0 "Y[N+1];" "6#&%\"YG6#,&%\"NG\"
\"\"F(F(" }{TEXT -1 49 " that we got from the boundary condition to yi
eld" }}{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 define
s " }{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 values at the 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(ode3,\{x=X(k)
,y(x)=y[k],\n diff(y(x),x)=Yp(k),diff(y(x),x$2)=Ypp(k)\} );\n
od;" }}}{PARA 0 "" 0 "" {TEXT -1 160 "We realize this system of N+1 eq
uations 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 p
oints." }}{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(lh
s(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 co
mes from the boundary condition 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 write this as a matrix. We list out the ro
ws 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 construc
t 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 u
nknown 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 61 "V:=Vector([rhs(bc3[1]),seq(rhs(eq[p
]),p=1..N-1),rhs(eq[N])]);" }}}{PARA 0 "" 0 "" {TEXT -1 48 "Finally, h
ere 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 sol
ve 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 "resul
t:=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 solut
ion." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "sol3:=combine(dsolve
(\{ode3,bc3\},y(x)));" }}}{PARA 0 "" 0 "" {TEXT -1 64 "Next, we superi
mpose the exact solution with this approximation." }}{EXCHG {PARA 0 ">
" 0 "" {MPLTEXT 1 0 115 "plot([rhs(sol3),result],x=-1..1,\n color=[
BLACK,GREEN], thickness=[1,3],\n legend=[`exact`,`finite difference`
]);" }}}{PARA 0 "" 0 "" {TEXT -1 81 "If you don't think this fit is pr
etty 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 methods, or discret
e replacement methods, for boundary value problems in differential equ
ations was made to suggest the techniques for partial differential equ
ations. We examine these next." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 50 "EMAIL: herod@math.gatech.edu or jhero
d@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" 88 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 1 1 2
33 1 1 }