{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 18 0 0 0 0 0 1 1 0 0 0 0 0 0 1 }{CSTYLE "" -1
257 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 258 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 259 "" 0 1 0 0 255 1 0 0 0 0 0
0 0 0 0 1 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
262 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "corrected" -1 263 ""
0 0 255 255 0 0 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 264 "" 0 1 0 0 0 0
0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 265 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0
0 1 }{CSTYLE "" -1 266 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE ""
-1 267 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 268 "" 0 1 0
0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 269 "" 0 1 0 0 0 0 1 0 0 0 0
0 0 0 0 1 }{CSTYLE "" -1 270 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 271 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
272 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "red" -1 273 "red" 0
0 255 0 0 1 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 274 "" 1 18 0 0 0 0 0
1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 275 "" 1 18 0 0 0 0 0 0 0 0 0 0 0 0
0 1 }{CSTYLE "" -1 276 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE ""
-1 277 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 278 "" 0 1
0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 279 "" 0 1 0 0 255 1 0 0 0
0 0 0 0 0 0 1 }{CSTYLE "" -1 280 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 281 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
282 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 283 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 284 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 285 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 286 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
287 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 288 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 289 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 290 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 291 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
292 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 293 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 294 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 295 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 296 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
297 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 298 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 299 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 300 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 301 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
302 "" 0 1 0 0 0 0 1 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 303 "" 0 1 0 0
0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 304 "" 0 1 0 0 0 0 0 1 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 305 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 306 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
307 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 308 "" 0 1 0 0
0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 309 "" 0 1 0 0 0 0 0 1 0 0 0 0
0 0 0 1 }{CSTYLE "" -1 310 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }
{CSTYLE "" -1 311 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1
312 "" 0 1 0 0 0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 313 "" 0 1 0 0
0 0 0 1 0 0 0 0 0 0 0 1 }{CSTYLE "" -1 314 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 0 }{CSTYLE "" -1 315 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 316 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1
317 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 318 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 319 "" 0 1 0 0 0 0 1 0 0 0 0 0
0 0 0 0 }{CSTYLE "" -1 320 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }
{CSTYLE "" -1 321 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1
322 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 323 "" 0 1 0 0
0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 324 "" 0 1 0 0 0 0 0 1 0 0 0 0
0 0 0 0 }{CSTYLE "" -1 325 "" 0 1 0 0 0 0 0 1 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 "Heading 2" -1 4 1 {CSTYLE "" -1 -1 "Times
" 1 14 0 0 0 1 2 1 2 2 2 2 1 1 1 1 }1 1 0 0 8 2 1 0 1 0 2 2 0 1 }
{PSTYLE "Heading 3" -1 5 1 {CSTYLE "" -1 -1 "Times" 1 12 0 0 0 1 1 1
2 2 2 2 1 1 1 1 }1 1 0 0 0 0 1 0 1 0 2 2 0 1 }{PSTYLE "Bullet Item"
-1 15 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 3 3 1 0 1 0 2 2 15 2 }{PSTYLE "Title" -1 18 1 {CSTYLE "" -1 -1 "
Times" 1 18 0 0 0 1 2 1 1 2 2 2 1 1 1 1 }3 1 0 0 12 12 1 0 1 0 2 2 19
1 }{PSTYLE "R3 Font 0" -1 256 1 {CSTYLE "" -1 -1 "Courier" 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 "R3 Font 2
" -1 257 1 {CSTYLE "" -1 -1 "Courier" 1 12 0 128 128 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 }}
{SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT 256 47 "Mathematical Models: Exa
mples for ODEs and PDEs" }{TEXT 274 0 "" }{TEXT 275 0 "" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT
-1 49 "Ordinary differential equations: population Model" }}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 187 "We want to model the growth of a populat
ion of individuals (e.g. rabbits). In a simple model, the growth of th
e population at time t would be proportional to the number of living i
ndividu" }{TEXT 263 3 "als" }{TEXT -1 354 " at that time (i.e. constan
t birth rate, model of Maltus). However, in this simple model the popu
lation would grow without bound. We want to consider a slightly better
model, in which we introduce a 'death rate' increasing with the numbe
r of living individuals (cf. slide 'Model Refinement 2' from lesson 4)
. In this model our population will be bounded." }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 8 "restart:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "with(DE
tools):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 32 "We assume a constant birth rate " }{TEXT
263 1 "b" }{TEXT -1 26 " (per living individuals)." }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 9 "b:=t->b0;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 23 "We
want the death rate " }{TEXT 263 1 "d" }{TEXT -1 42 " to increase lin
early with the population:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 17 "d:=t-
>d0+d1*p(t);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 82 "We get the follow
ing ordinary differential equation (ODE) for the population p(t):" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "dgl:=diff(p(t),t)=(b(t)-d(t))*p(t);
" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 24 "We choose the constants " }
{TEXT 263 2 "b0" }{TEXT -1 27 ", d0, d1 and p0 in our ODE " }{TEXT
258 3 "dgl" }{TEXT -1 1 ":" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "smpld
ata := \{b0=0.07, d0=0.04, d1=0.0005, p0=6\};" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 25 "dgl1:=subs(smpldata,dgl);" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 41 "This ODE is sometimes written in the form" }}{PARA 0 ""
0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 24 "p'(t) = a p(t) - b p
(t)\262" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1
18 "and is called the " }{TEXT 257 17 "Logistic Equation" }{TEXT -1 1
"." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT
-1 113 "The ODE is very simple (as it can be separated). Maple can sol
ve it directly. Note that we use the initial value " }{TEXT 265 7 "p(0
)=p0" }{TEXT -1 29 " to obtain a unique solution." }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 34 "sol := dsolve(\{dgl,p(0)=p0\},p(t));" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "pex := unapply(subs(sol,p(t)),t);"
}}{PARA 0 "" 0 "" {TEXT 263 50 "(unapply converts our solution to a fu
nction of t)" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 32 "We can calculate \+
the limit t -> " }{XPPEDIT 18 0 "infinity" "6#%)infinityG" }}{PARA 0 "
> " 0 "" {MPLTEXT 1 0 41 "limit (subs(smpldata,pex(t)),t=infinity);" }
}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 131 "Before we proceed, we want to d
iscuss some basics about ordinary differential equations. In general a
first order ODE has the form:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}
{PARA 0 "" 0 "" {TEXT -1 4 "Let " }{TEXT 260 9 "f: G -> R" }{TEXT -1
92 " be a real continuous function on a subset G of the real plane R
\262. We consider the equation" }}{PARA 0 "" 0 "" {TEXT 261 12 "y' = \+
f (t,y)" }{TEXT 264 1 " " }{TEXT -1 6 " " }}{PARA 0 "" 0 ""
{TEXT -1 74 "and are interested in differentiable functions y(t) for w
hich the equation" }}{PARA 0 "" 0 "" {TEXT 262 17 "y'(t) = f(t,y(t))"
}{TEXT -1 1 " " }}{PARA 0 "" 0 "" {TEXT -1 18 "holds for every t." }}
{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 68 "A geometr
ic interpretation of an ODE is the following. The equation " }{TEXT
267 9 "y'=f(t,y)" }{TEXT -1 43 " defines a direction field. In every p
oint " }{TEXT 268 5 "(t,y)" }{TEXT -1 128 " a tangent for y(t) is give
n. Every curve y(t) that fits into the direction field is a solution o
f the ODE. You can use Maple's " }{TEXT 259 9 "DEplot( )" }{TEXT -1
89 " command to plot the direction field of an ODE. Note that all curv
es seem to converge to " }{TEXT 266 4 "p=60" }{TEXT -1 50 ". This corr
esponds to the limit calculation above." }}{PARA 0 "> " 0 "" {MPLTEXT
1 0 36 "DEplot(dgl1,p(t),t=0..250,p=0..100);" }}}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 86 "If you look at the plot above, it is clear that there i
s an infinite number of curves " }{TEXT 263 4 "p(t)" }{TEXT -1 100 " t
hat fit into the direction field. But we will obtain a unique solution
if we fix an initial value " }{TEXT 269 8 "p(t0)=p0" }{TEXT -1 12 ". \+
We choose " }{TEXT 270 6 "p(0)=6" }{TEXT -1 100 " and plot the directi
on field again, this time with the unique solution of the inital value
problem:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "DEplot(dgl1,p(
t),t=0..250,\{[0,6]\},p=0..100);" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "Given an init
ial value problem " }{TEXT 271 9 "y'=f(x,y)" }{TEXT -1 2 ", " }{TEXT
272 9 "y(t0)=y0," }{TEXT -1 171 " two important questions arise: The q
uestion of existence of solutions and if so the question of uniqueness
. We don't want to go into details but only state the following:" }}
{PARA 0 "" 0 "" {TEXT -1 236 "A sufficient condition that there is a s
olution for the initial value problem is, that the function f is 'Lips
chitz-continuous'. This is known as the 'Picard-Lindel\366f-Theorem'.
Moreover, under this circumstances the solution is unique." }}{PARA
0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 102 "For our initia
l value problem above the 'Lipschitz-condition' is satisfied and the s
olution is unique." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}
{SECT 1 {PARA 3 "" 0 "" {TEXT -1 33 "Exercise: Verhulst, Predator-Prey
" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}{PARA 0 "> "
0 "" {MPLTEXT 1 0 12 "with(plots):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
14 "with(DEtools):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{SECT 1 {PARA 4 "" 0 "" {TEXT -1 17 "Model of Verhulst" }}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 12 "with(plots):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "wi
th(DEtools):" }}{PARA 0 "" 0 "" {TEXT -1 50 "Use Maple to solve the si
mpler model of Verhulst: " }{XPPEDIT 18 0 "diff(p(t),t) = gamma[0]-del
ta[0]-(gamma[1]+delta[1])*p(t);" "6#/-%%diffG6$-%\"pG6#%\"tGF*,(&%&gam
maG6#\"\"!\"\"\"&%&deltaG6#F/!\"\"*&,&&F-6#F0F0&F26#F0F0F0-F(6#F*F0F4
" }{TEXT -1 18 ", with parameters " }{XPPEDIT 18 0 "gamma[0]-delta[0] \+
= 6;" "6#/,&&%&gammaG6#\"\"!\"\"\"&%&deltaG6#F(!\"\"\"\"'" }{TEXT -1
5 " and " }{XPPEDIT 18 0 "gamma[1]+delta[1] = 1/10;" "6#/,&&%&gammaG6#
\"\"\"F(&%&deltaG6#F(F(*&F(F(\"#5!\"\"" }{TEXT -1 1 "." }}{PARA 0 ""
0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 61 "First, store
the differential equation in a variable called " }{TEXT 0 4 "verh" }
{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA
0 "" 0 "" {TEXT -1 142 "Next, use DEplot to draw the direction field o
f this ODE (reasonable ranges are t=0..100 and p=0.120) with a solutio
n curve starting at p(0)=6" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 84 "For this simple ODE, we can easily
compute exact solutions - do this (using Maple's " }{TEXT 0 6 "dsolve
" }{TEXT -1 2 "):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT
-1 25 "Predator-Prey simulation " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1
146 "Next, we consider a population model with two species, a predato
r (lions, e.g.) and a prey (gazelles or whatever lions may eat). Now t
he numbers " }{XPPEDIT 18 0 "p(t);" "6#-%\"pG6#%\"tG" }{TEXT -1 18 " o
f predators and " }{XPPEDIT 18 0 "q(t);" "6#-%\"qG6#%\"tG" }{TEXT -1
108 " of prey relate to the respective growth rate. The Lotka-Volterra
model assumes a relationship like follows:" }}{PARA 0 "" 0 ""
{XPPEDIT 18 0 "diff(p(t),t) = (b[1]+c[1]*q)*p;" "6#/-%%diffG6$-%\"pG6#
%\"tGF**&,&&%\"bG6#\"\"\"F0*&&%\"cG6#F0F0%\"qGF0F0F0F(F0" }{TEXT -1 4
" and" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }{XPPEDIT 18 0 "diff(q(t),t) = \+
(b[2]+c[2]*p)*q;" "6#/-%%diffG6$-%\"qG6#%\"tGF**&,&&%\"bG6#\"\"#\"\"\"
*&&%\"cG6#F0F1%\"pGF1F1F1F(F1" }}{PARA 0 "" 0 "" {TEXT -1 5 "with " }
{XPPEDIT 18 0 "0 < c1;" "6#2\"\"!%#c1G" }{TEXT -1 36 " (predators need
something to eat), " }{XPPEDIT 18 0 "c2 < 0;" "6#2%#c2G\"\"!" }{TEXT
-1 33 " (more predators eat more prey), " }{XPPEDIT 18 0 "b1 < 0;" "6#
2%#b1G\"\"!" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "0 < b2;" "6#2\"\"!%#b
2G" }{TEXT -1 1 "." }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restar
t;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}{PARA 0 "> "
0 "" {MPLTEXT 1 0 14 "with(DEtools):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT
-1 47 "Store the two equations above in the variables " }{TEXT 0 2 "L1
" }{TEXT -1 5 " and " }{TEXT 0 2 "L2" }{TEXT -1 15 ", respectively:" }
}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 13 "Define a set " }{TEXT 0 3
"sys" }{TEXT -1 36 " of the two equations, substituting " }{TEXT 0 7 "
b1=-1/2" }{TEXT -1 2 ", " }{TEXT 0 6 "b2=1/5" }{TEXT -1 2 ", " }{TEXT
0 8 "c1=1/200" }{TEXT -1 2 ", " }{TEXT 0 8 "c2=-1/50" }{TEXT -1 118 " \+
(of course, you can choose parameters of your own, if you like, but I \+
recommend starting with the values given here):" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 52 "Draw a plot o
f the directional field for the system " }{TEXT 0 3 "sys" }{TEXT -1 7
" using " }{TEXT 0 6 "DEplot" }{TEXT -1 65 " (Hint: you have to specif
y a range for the independent variable " }{TEXT 320 1 "t" }{TEXT -1
74 ", even if this range does not affect the plot as the axes now are \+
labeled " }{TEXT 319 1 "p" }{TEXT -1 5 " and " }{TEXT 321 1 "q" }
{TEXT -1 57 ". For the coefficients from above, reasonable ranges for \+
" }{TEXT 322 1 "p" }{TEXT -1 5 " and " }{TEXT 323 1 "q" }{TEXT -1 24 "
are p=0..30, q=0..200):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 96 "Maple is not able to give us the e
xact solution of this set of ODEs but it (i.e. its procedure " }
{TEXT 0 6 "DEplot" }{TEXT -1 122 ") can compute solution curves for gi
ven initial values numerically (we will learn more about this in the n
ext worksheet). " }}{PARA 15 "" 0 "" {TEXT -1 284 "You have to specify
a set of initial conditions, say p(0)=6, q(0)=50. In the call to DEpl
ot, they are specified between the range for t and the ranges for p an
d q as a list containing a list of the two conditions (which simply me
ans: with a double pair of brackets [...] around them)." }}{PARA 15 "
" 0 "" {TEXT -1 85 "Now the range for t affects the plot (in which way
?). A reasonable range is t=0..100." }}{PARA 15 "" 0 "" {TEXT -1 117 "
You should specify a stepsize for the numerical computation - a stepsi
ze of 0.2 will result in a fairly smooth curve." }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 121 "You can use \+
the additional keyword scene=[t,p] in DEplot to get a plot of the numb
er of predators versus the time axis t." }}{PARA 0 "" 0 "" {TEXT -1
146 "Do this (using the procedure display) to produce a plot containin
g the two curves 'number of predators' and 'number of prey' versus the
time axis." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 27 "So
lution: Model of Verhulst" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "
restart;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 14 "with(DEtools):" }}}{PARA 0 "" 0 "" {TEXT
-1 50 "Use Maple to solve the simpler model of Verhulst: " }{XPPEDIT
18 0 "diff(p(t),t) = gamma[0]-delta[0]-(gamma[1]+delta[1])*p(t);" "6#/
-%%diffG6$-%\"pG6#%\"tGF*,(&%&gammaG6#\"\"!\"\"\"&%&deltaG6#F/!\"\"*&,
&&F-6#F0F0&F26#F0F0F0-F(6#F*F0F4" }{TEXT -1 18 ", with parameters " }
{XPPEDIT 18 0 "gamma[0]-delta[0] = 6;" "6#/,&&%&gammaG6#\"\"!\"\"\"&%&
deltaG6#F(!\"\"\"\"'" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "gamma[1]+del
ta[1] = 1/10;" "6#/,&&%&gammaG6#\"\"\"F(&%&deltaG6#F(F(*&F(F(\"#5!\"\"
" }{TEXT -1 1 "." }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 61 "First, store t
he differential equation in a variable called " }{TEXT 0 4 "verh" }
{TEXT -1 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "verh := diff(p(t),
t) = 6 - 0.1*p(t);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 142 "Next, use \+
DEplot to draw the direction field of this ODE (reasonable ranges are \+
t=0..100 and p=0.120) with a solution curve starting at p(0)=6" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 48 "DEplot(verh, p(t), t=0..100, \{[0,6
]\}, p=0..120);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 84 "For this simpl
e ODE, we can easily compute exact solutions - do this (using Maple's \+
" }{TEXT 0 6 "dsolve" }{TEXT -1 2 "):" }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 19 "dsolve(verh, p(t));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0
"" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 23 "Solution: Predator-Prey" }
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 146 "Next, we consider a population \+
model with two species, a predator (lions, e.g.) and a prey (gazelles \+
or whatever lions may eat). Now the numbers " }{XPPEDIT 18 0 "p(t);" "
6#-%\"pG6#%\"tG" }{TEXT -1 18 " of predators and " }{XPPEDIT 18 0 "q(t
);" "6#-%\"qG6#%\"tG" }{TEXT -1 108 " of prey relate to the respective
growth rate. The Lotka-Volterra model assumes a relationship like fol
lows:" }}{PARA 0 "" 0 "" {XPPEDIT 18 0 "diff(p(t),t) = (b[1]+c[1]*q)*p
;" "6#/-%%diffG6$-%\"pG6#%\"tGF**&,&&%\"bG6#\"\"\"F0*&&%\"cG6#F0F0%\"q
GF0F0F0F(F0" }{TEXT -1 4 " and" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }
{XPPEDIT 18 0 "diff(q(t),t) = (b[2]+c[2]*p)*q;" "6#/-%%diffG6$-%\"qG6#
%\"tGF**&,&&%\"bG6#\"\"#\"\"\"*&&%\"cG6#F0F1%\"pGF1F1F1F(F1" }}{PARA
0 "" 0 "" {TEXT -1 5 "with " }{XPPEDIT 18 0 "0 < c1;" "6#2\"\"!%#c1G"
}{TEXT -1 36 " (predators need something to eat), " }{XPPEDIT 18 0 "c2
< 0;" "6#2%#c2G\"\"!" }{TEXT -1 33 " (more predators eat more prey), \+
" }{XPPEDIT 18 0 "b1 < 0;" "6#2%#b1G\"\"!" }{TEXT -1 5 " and " }
{XPPEDIT 18 0 "0 < b2;" "6#2\"\"!%#b2G" }{TEXT -1 1 "." }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 12 "with(plots):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "wi
th(DEtools):" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 47 "Store the two equ
ations above in the variables " }{TEXT 0 2 "L1" }{TEXT -1 5 " and " }
{TEXT 0 2 "L2" }{TEXT -1 15 ", respectively:" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 39 "L1 := diff(p(t),t) = (b1+c1*q(t))*p(t);" }}{PARA 0 ">
" 0 "" {MPLTEXT 1 0 39 "L2 := diff(q(t),t) = (b2+c2*p(t))*q(t);" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 13 "Define a set " }{TEXT 0 3 "sys" }
{TEXT -1 36 " of the two equations, substituting " }{TEXT 0 7 "b1=-1/2
" }{TEXT -1 2 ", " }{TEXT 0 6 "b2=1/5" }{TEXT -1 2 ", " }{TEXT 0 8 "c1
=1/200" }{TEXT -1 2 ", " }{TEXT 0 8 "c2=-1/50" }{TEXT -1 118 " (of cou
rse, you can choose parameters of your own, if you like, but I recomme
nd starting with the values given here):" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 56 "sys := subs(\{b1=-1/2,b2=1/5,c1=1/200,c2=-1/50\},\{L1
,L2\});" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 52 "Draw a plot of the dir
ectional field for the system " }{TEXT 0 3 "sys" }{TEXT -1 7 " using \+
" }{TEXT 0 6 "DEplot" }{TEXT -1 65 " (Hint: you have to specify a rang
e for the independent variable " }{TEXT 315 1 "t" }{TEXT -1 74 ", even
if this range does not affect the plot as the axes now are labeled "
}{TEXT 314 1 "p" }{TEXT -1 5 " and " }{TEXT 316 1 "q" }{TEXT -1 57 ". \+
For the coefficients from above, reasonable ranges for " }{TEXT 317 1
"p" }{TEXT -1 5 " and " }{TEXT 318 1 "q" }{TEXT -1 24 " are p=0..30, q
=0..200):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 43 "DEplot(sys,[p,q],t=0..
50,p=0..30,q=0..200);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 96 "Maple is
not able to give us the exact solution of this set of ODEs but it (i
.e. its procedure " }{TEXT 0 6 "DEplot" }{TEXT -1 122 ") can compute s
olution curves for given initial values numerically (we will learn mor
e about this in the next worksheet). " }}{PARA 15 "" 0 "" {TEXT -1
284 "You have to specify a set of initial conditions, say p(0)=6, q(0)
=50. In the call to DEplot, they are specified between the range for t
and the ranges for p and q as a list containing a list of the two con
ditions (which simply means: with a double pair of brackets [...] arou
nd them)." }}{PARA 15 "" 0 "" {TEXT -1 85 "Now the range for t affects
the plot (in which way?). A reasonable range is t=0..100." }}{PARA
15 "" 0 "" {TEXT -1 117 "You should specify a stepsize for the numeric
al computation - a stepsize of 0.2 will result in a fairly smooth curv
e." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 75 "DEplot(sys,[p,q],t=0..50,[[p(
0)=6,q(0)=50]],p=0..30,q=0..200,stepsize=0.2);" }}}{EXCHG {PARA 0 ""
0 "" {TEXT -1 121 "You can use the additional keyword scene=[t,p] in D
Eplot to get a plot of the number of predators versus the time axis t.
" }}{PARA 0 "" 0 "" {TEXT -1 146 "Do this (using the procedure display
) to produce a plot containing the two curves 'number of predators' an
d 'number of prey' versus the time axis." }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 93 "p1 := DEplot(sys,[p,q],t=0..100,[[p(0)=6,q(0)=50]],st
epsize=0.2, scene=[t,p], linecolor=red):" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 93 "p2 := DEplot(sys,[p,q],t=0..100,[[p(0)=6,q(0)=50]],st
epsize=0.2, scene=[t,q],linecolor=blue):" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 17 "display([p1,p2]);" }}}{EXCHG {PARA 0 "> " 0 ""
{MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}
{SECT 1 {PARA 3 "" 0 "" {TEXT -1 47 "Partial differential equations: h
eat conduction" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 36 "The one dimensi
onal heat conduction " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 295 "We want \+
to consider the problem of heat conducting in a medium (without curren
ts or radiation) in the one dimensional case. We are given a wire (of \+
infinite length) which has a given distribution of temperature at time
t=0. The question is how the heat is conducted through the body of th
e wire." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1
46 "The situation can be modeled by the so called " }{TEXT 276 13 "hea
t equation" }{TEXT -1 144 ". It is an example of a partial differentia
l equation (PDE). If u(x,t) is the temperature at position x and time \+
t the equation looks like this:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}}{EXCHG {PARA 0 "
> " 0 "" {MPLTEXT 1 0 42 "heat := diff(u(x,t),t)=k*diff(u(x,t),x,x);"
}}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1
295 "The constant k depends on the heat conductivity of the medium. Th
e heat equation shows what a PDE is in general: We are given some kind
of relationship between a function u(x,t) and/or its partial derivati
ves. In the one dimensional case this can be partial derivatives in ti
me t or in space x. " }}{PARA 0 "" 0 "" {TEXT -1 172 "The problem is t
o find functions u(x,t) which satisfy the given relationship and somet
imes satisfy additional conditions (initital value conditions or bound
ary conditions)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 ""
{TEXT -1 270 "PDEs like the heat equation (but usually much more compl
ex) very often arise in modeling problems in science, engineering and \+
financial engineering. Sometimes we are able to find exact solutions u
sing analytic methods. Maple is also capable of solving some PDEs via \+
the " }{TEXT 277 11 "pdesolve( )" }{TEXT -1 16 " command (Click " }
{HYPERLNK 17 "here" 2 "pdesolve" "" }{TEXT -1 199 " for details of thi
s command). However, often a PDE cannot be solved analytically and our
only chance is a numerical (and therefore approximate) solution (we w
ill discuss that later in this lecture)." }}{PARA 0 "" 0 "" {TEXT -1
0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 132 "In this exercise we want \+
to solve the heat equation exactly. Older Versions of Maple cannot sol
ve it automatically, but Maple 6 can:" }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 22 "pdesolve(heat,u(x,t));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 0 "" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 22 "Solution by separation
" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 69 "A typical technique to solve t
his type of equations without Maple is " }{TEXT 278 13 "by separation
" }{TEXT -1 79 ". This means we assume that our solution u(x,t) is of \+
the form u(x,t)=X(x)T(t)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "eq := \+
subs(u(x,t)=X(x)*T(t),heat);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 86 "A
s X(x) is independent of t and T(t) of x the partial derivatives beco
me much easier." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "sep := eq/X(x)/T
(t);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 290 "We now have \"separated
\" the variables: each side of the equation depends only on one variab
le. The only way this equation can be satisfied is if both sides are c
onstant. We only consider the case of negative values (this is the cas
e that is of physical interest) and denote the constant by " }
{XPPEDIT 18 0 "-omega^2;" "6#,$*$%&omegaG\"\"#!\"\"" }{TEXT -1 53 ". S
o we have for the left hand side of the equation: " }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 19 "lhs(sep)= -omega^2;" }{TEXT -1 0 "" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 138 "This is an ordinary differential equatio
n (as T(t) is a function of one variable). This equation can easily be
solved. We can use Maple's " }{TEXT 279 10 "dsolve( ) " }{TEXT -1 33
"command and obtain the following:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
23 "Tsol := dsolve(%,T(t));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 61 "Th
e same can be done for the right hand side of the equation:" }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 19 "rhs(sep)= -omega^2;" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 23 "Xsol := dsolve(%,X(x));" }{TEXT -1 0 "" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 40 "We now have solutions X(x) and T(t). In \+
" }{XPPEDIT 18 0 "Tsol;" "6#%%TsolG" }{TEXT -1 43 " we have the consta
nt (due to integration) " }{TEXT 283 3 "_C1" }{TEXT -1 8 " and in " }
{XPPEDIT 18 0 "Xsol;" "6#%%XsolG" }{TEXT -1 15 " the constants " }
{TEXT 284 3 "_C1" }{TEXT -1 6 " and " }{TEXT 285 6 "_C2. " }{TEXT
-1 15 "The constants _" }{TEXT 280 2 "C1" }{TEXT -1 4 " in " }
{XPPEDIT 18 0 "Tsol" "6#%%TsolG" }{TEXT -1 6 " and _" }{TEXT 281 3 "C1
" }{TEXT -1 3 "in " }{XPPEDIT 18 0 "Xsol" "6#%%XsolG" }{TEXT 282 1 " \+
" }{TEXT -1 78 "are not the same of course, so we need to rename them \+
(contrary, the constant " }{XPPEDIT 18 0 "omega;" "6#%&omegaG" }{TEXT
286 2 " " }{TEXT 287 3 "is " }{TEXT -1 56 " really the same constant \+
due to our separation ansatz):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "X
sol := subs(_C2=_C3,_C1=_C2,%);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1
50 "The complete solution of our problem therefore is " }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 26 "usol:=rhs(Tsol)*rhs(Xsol);" }}}{EXCHG {PARA 0
"" 0 "" {TEXT -1 14 "The constants " }{XPPEDIT 18 0 "omega;" "6#%&omeg
aG" }{TEXT 288 15 ", _C1, _C2, _C3" }{TEXT -1 339 " can be chosen arbi
trarily. Usually we are given some initial value conditions e.g. the v
alues u(x,t=0) which would in our case be the temperature distribution
at t=0. Using this additional information we can calculate the consta
nts and obtain a unique solution. In this exercise we just choose the \+
constants and the heat conductivity k by" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 57 "usol_par := subs(_C1=1/2,_C2=1,_C3=1/2,omega=1,k=4,us
ol);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 29 "Now we can plot the solut
ion." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "plot3d(usol_par,x=-2*Pi..2*
Pi,t=0..10,axes=boxed);" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 185 "Maple offers a nice feature: We can animate our result t
o see how the heat is conducting. You might have to right-click on the
graphic and use the menu option 'animation-play' to start." }}{PARA
0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 61 "animate(usol_par,x=-2*Pi..2*Pi,t=0..10,frames=100,col
or=red);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 32 "Finally we test if ou
r solution " }{TEXT 289 5 "usol " }{TEXT -1 18 "satisfies the PDE:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 23 "subs(u(x,t)=usol,heat);" }}}{EXCHG
{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "simplify(%);" }}}{EXCHG {PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 14 "Four
ier method" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "In the previous sect
ion, we computed solutions for the heat conduction equation of the typ
e " }{XPPEDIT 18 0 "_C1*exp(-omega^2*t)*(_C2*cos(omega*x/sqrt(k))+_C3*
sin(omega*x/sqrt(k)));" "6#*(%$_C1G\"\"\"-%$expG6#,$*&%&omegaG\"\"#%\"
tGF%!\"\"F%,&*&%$_C2GF%-%$cosG6#*(F+F%%\"xGF%-%%sqrtG6#%\"kGF.F%F%*&%$
_C3GF%-%$sinG6#*(F+F%F6F%-F86#F:F.F%F%F%" }}{PARA 0 "" 0 "" {TEXT -1
41 "that resulted from our separation ansatz " }{XPPEDIT 18 0 "u(x,t) \+
= X(x)*T(t);" "6#/-%\"uG6$%\"xG%\"tG*&-%\"XG6#F'\"\"\"-%\"TG6#F(F-" }
{TEXT -1 2 ". " }}{PARA 0 "" 0 "" {TEXT -1 197 "Now, we will extend th
is to more general inital temerature distributions by the observation \+
that the sum of solutions of this form does also solve the equation (s
ince it is linear and homogeneous)." }}}{EXCHG {PARA 0 "" 0 "" {TEXT
-1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 25 "We make use of this by a " }
{TEXT 290 23 "Fourier Series approach" }{TEXT -1 79 ". The idea behind
it is that we can represent under fairly general assumptions " }
{XPPEDIT 18 0 "2*Pi;" "6#*&\"\"#\"\"\"%#PiGF%" }{TEXT -1 31 "-periodic
ally functions by its " }{TEXT 291 17 "Fourier Series. " }{TEXT -1
85 "Therefore we have for a Fourier transformation with respect to the
spatial direction " }{XPPEDIT 18 0 "x;" "6#%\"xG" }{TEXT -1 8 " (with
" }{XPPEDIT 18 0 "a(omega,t);" "6#-%\"aG6$%&omegaG%\"tG" }{TEXT -1
46 " as the time dependent fourier coefficients):" }}{PARA 0 "" 0 ""
{XPPEDIT 18 0 "u(x,t) = 1/sqrt(pi);" "6#/-%\"uG6$%\"xG%\"tG*&\"\"\"F*-
%%sqrtG6#%#piG!\"\"" }{TEXT 292 2 " " }{XPPEDIT 18 0 "sum(a(omega,t)*
exp(I*omega*x),omega = 0 .. infinity) = 1/sqrt(Pi);" "6#/-%$sumG6$*&-%
\"aG6$%&omegaG%\"tG\"\"\"-%$expG6#*(%\"IGF-F+F-%\"xGF-F-/F+;\"\"!%)inf
inityG*&F-F--%%sqrtG6#%#PiG!\"\"" }{XPPEDIT 18 0 "sum(a(omega,t)*(cos(
omega*t)+I*sin(omega*t)),omega = 0 .. infinity);" "6#-%$sumG6$*&-%\"aG
6$%&omegaG%\"tG\"\"\",&-%$cosG6#*&F*F,F+F,F,*&%\"IGF,-%$sinG6#*&F*F,F+
F,F,F,F,/F*;\"\"!%)infinityG" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1
0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 61 "If we substitute this fo
rm for u(x,t) into the heat equation " }{XPPEDIT 18 0 "diff(u(x,t),t) \+
= k*diff(u(x,t),`$`(x,2));" "6#/-%%diffG6$-%\"uG6$%\"xG%\"tGF+*&%\"kG
\"\"\"-F%6$-F(6$F*F+-%\"$G6$F*\"\"#F." }{TEXT -1 80 ", we see that eac
h term of the sum must satisfy the PDE and therefore, we obtain" }}
{PARA 0 "" 0 "" {XPPEDIT 18 0 "a(omega,t) = a(omega)*exp(-omega^2*t);
" "6#/-%\"aG6$%&omegaG%\"tG*&-F%6#F'\"\"\"-%$expG6#,$*&F'\"\"#F(F,!\"
\"F," }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 55 "w
hich lead again the well-known solutions, if we write " }{XPPEDIT 18
0 "exp(iwx)" "6#-%$expG6#%$iwxG" }{TEXT -1 13 " in terms of " }{TEXT
294 3 "sin" }{TEXT -1 5 " and " }{TEXT 293 6 "cos. " }{TEXT -1 121 "F
or simplicity, in the following we will allow only odd functions (whi
ch means we can forget the imaginary part and all " }{TEXT 295 3 "cos
" }{TEXT -1 17 " terms), and get:" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }
{XPPEDIT 18 0 "u(x,t) = 1/sqrt(pi);" "6#/-%\"uG6$%\"xG%\"tG*&\"\"\"F*-
%%sqrtG6#%#piG!\"\"" }{TEXT 296 2 " " }{XPPEDIT 18 0 "sum(a(omega)*ex
p(-omega^2*t)*sin(omega*x),omega = 1 .. infinity);" "6#-%$sumG6$*(-%\"
aG6#%&omegaG\"\"\"-%$expG6#,$*&F*\"\"#%\"tGF+!\"\"F+-%$sinG6#*&F*F+%\"
xGF+F+/F*;F+%)infinityG" }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 105 "So we can compute the time evolution of \+
arbitrary initial temperature distributions u(x,0), provided that" }}
{PARA 15 "" 0 "" {XPPEDIT 18 0 "u(x,0);" "6#-%\"uG6$%\"xG\"\"!" }
{TEXT -1 6 " is a " }{XPPEDIT 18 0 "2*Pi;" "6#*&\"\"#\"\"\"%#PiGF%" }
{TEXT -1 22 "-periodical function: " }{XPPEDIT 18 0 "u(x+2*Pi,0) = u(x
,0);" "6#/-%\"uG6$,&%\"xG\"\"\"*&\"\"#F)%#PiGF)F)\"\"!-F%6$F(F-" }
{TEXT -1 1 "," }}{PARA 15 "" 0 "" {XPPEDIT 18 0 "u(x,0);" "6#-%\"uG6$%
\"xG\"\"!" }{TEXT -1 41 " is an odd function (with respect to x): " }
{XPPEDIT 18 0 "u(-x) = -u(x);" "6#/-%\"uG6#,$%\"xG!\"\",$-F%6#F(F)" }
{TEXT -1 19 ", which means that " }{XPPEDIT 18 0 "u(x,0);" "6#-%\"uG6$
%\"xG\"\"!" }{TEXT -1 47 " can be expressed as a series of sin-functio
ns," }}{PARA 15 "" 0 "" {TEXT -1 80 "and the Fourier series does conve
rge to u (which is only a very weak condition)." }}{PARA 0 "" 0 ""
{TEXT -1 182 "So expanding our initial temperature distribution u(x,t=
0) in terms of the functions above is reasonable. Of course, Maple can
not store the infinite number of arbitrary coefficients " }{XPPEDIT
18 0 "a(omega),omega = 1 .. infinity;" "6$-%\"aG6#%&omegaG/F&;\"\"\"%)
infinityG" }{TEXT -1 28 ", so we only use the first " }{TEXT 297 5 "n
max " }{TEXT -1 68 "terms of our fourier expansion. First let Maple ca
lculate the terms " }{XPPEDIT 18 0 "exp(-omega^2*t)*sin(omega*x);" "6#
*&-%$expG6#,$*&%&omegaG\"\"#%\"tG\"\"\"!\"\"F,-%$sinG6#*&F)F,%\"xGF,F,
" }{TEXT -1 47 " for fixed t=0 (giving a list of expressions " }
{TEXT 298 3 "fkn" }{TEXT -1 22 ") and for variable t (" }{TEXT 299 4 "
fknt" }{TEXT -1 2 ")." }}{PARA 0 "" 0 "" {TEXT -1 25 "Remember that th
e factor " }{XPPEDIT 18 0 "1/sqrt(Pi);" "6#*&\"\"\"F$-%%sqrtG6#%#PiG!
\"\"" }{TEXT -1 43 " is necessary for the functions to be ortho" }
{TEXT 300 6 "normal" }{TEXT -1 21 " with respect to the " }{XPPEDIT
18 0 "L[2];" "6#&%\"LG6#\"\"#" }{TEXT -1 20 "-scalar product (in " }
{XPPEDIT 18 0 "0 .. 2*Pi;" "6#;\"\"!*&\"\"#\"\"\"%#PiGF'" }{TEXT -1 2
")." }}}{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 10 "nmax := 5;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
68 "fknt := [seq(sin(omega*x)*exp(-omega^2*t)/sqrt(Pi), omega=1..nmax)
];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "fkn := simplify(subs(t=0,fknt
));" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 51 "Let's plot \+
these functions to see what we have got:" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 20 "plot(fkn,x=0..2*Pi);" }{TEXT -1 0 "" }}}{EXCHG {PARA
0 "> " 0 "" {MPLTEXT 1 0 53 "animate(\{seq(fknt[i],i=1..nmax)\}, x=0..
2*Pi, t=0..1);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 98 "As already mentioned, these functi
ons must satisfy the heat equation. Maple can show this quickly:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "seq(diff(i,x,x)-diff(i,t),i=fknt);
" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}
{SECT 1 {PARA 5 "" 0 "" {TEXT -1 67 "Expansion of an initial temperatu
re distribution: A simple example " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1
150 "Now we've everything we need to start. Let's choose as an initial
temperature distribution u(x,t=0) a simple combination of three fouri
er components. " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "w := [se
q(0,i=1..nops(fknt))];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "w[1] := -
1; w[2]:=-4; w[4]:=2;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 42 "ww := [seq
(w[i]*fknt[i],i=1..nops(fknt))];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1
49 "Now we've the fourier coefficients in the vector " }{TEXT 309 1 "w
" }{TEXT -1 38 " and the fourier components in vector " }{TEXT 308 2 "
ww" }{TEXT -1 106 ". We now only have to sum over all fourier componen
ts to obtain our initial temperature distribution h(x)." }}{PARA 0 "> \+
" 0 "" {MPLTEXT 1 0 32 "h := add(ww[j],j=1..nops(fknt));" }}{PARA 0 "
" 0 "" {TEXT -1 15 "A plot showing " }{TEXT 310 4 "h(x)" }{TEXT -1 28
" and its fourier components:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "pl
ot([subs(t=0,h),seq(subs(t=0,ww[i]),i=1..nops(fknt))] ,x=0..2*Pi);" }}
}{EXCHG {PARA 0 "" 0 "" {TEXT -1 49 "Now we can look how the heat of o
ur distribution " }{TEXT 311 4 "h(x)" }{TEXT -1 16 " is conducting. "
}}{PARA 0 "" 0 "" {TEXT -1 25 "First without components:" }{MPLTEXT 1
0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 62 "with(plots):animate(\{h\},
x=0..2*Pi,t=0..2,frames=32,color=red);" }{TEXT -1 0 "" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 41 "And together with its fourier components:
" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 69 "with(plots):animate(\{h,op(ww)
\},x=0..2*Pi,t=0..2,frames=32,color=red);" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 100 "You can see that the high frequency components decay muc
h faster than the low frequency components. " }{MPLTEXT 1 0 0 "" }}}
{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}}{SECT 1 {PARA 3 "" 0
"" {TEXT -1 67 "Exercise: The sawtooth function as initial temperature
distribution" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 101 "Now it is your t
urn. We want to choose as an initial temperature distribution u(x,t=0)
the so called " }{TEXT 306 18 "sawtooth function:" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 8 "restart;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 15 "g := \+
x -> x-Pi;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "plot(g(x mod 2*Pi) ,x
=-2*Pi..6*Pi);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 81 "First we redefi
ne the fourier functions (because we did a 'restart' to make sure)" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "nmax := 5;" }}{PARA 0 "> " 0 ""
{MPLTEXT 1 0 68 "fknt := [seq(sin(omega*x)*exp(-omega^2*t)/sqrt(Pi), o
mega=1..nmax)];" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "fkn := simplify(
subs(t=0,fknt));" }{TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1
10 "Note that " }{TEXT 0 6 "fkn[i]" }{TEXT -1 30 " is not a Maple func
tion, but " }{TEXT 0 17 "unapply(fkn[i],x)" }{TEXT -1 8 " is one." }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT
-1 14 "Your Workspace" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "We now ne
ed to find the Fourier representation of g. More precise, we have to c
alculate the " }{TEXT 307 21 "fourier coefficients " }{XPPEDIT 18 0 "i
nt(g(x)*fkn[i](x),x = 0 .. 2*Pi);" "6#-%$intG6$*&-%\"gG6#%\"xG\"\"\"-&
%$fknG6#%\"iG6#F*F+/F*;\"\"!*&\"\"#F+%#PiGF+" }{TEXT -1 37 " with our
functions fkn[i=1..nmax]. " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 62 "Le
t Maple calculate them and store the result in a list named " }{TEXT
301 2 "c!" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }{TEXT -1 0 "" }}}
{EXCHG {PARA 0 "" 0 "" {TEXT -1 23 "Now build a list named " }{TEXT
303 3 "cc " }{TEXT -1 30 "which contains the functions " }{TEXT 324
4 "fknt" }{TEXT -1 58 " multiplied by their coefficients from our four
ier series:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "
" 0 "" {TEXT -1 120 "We now can sum up all these terms and obtain an a
pproximation of the sawtooth function. We want to name this expression
" }{TEXT 304 2 "gg" }{TEXT -1 2 ". " }}{PARA 0 "> " 0 "" {MPLTEXT 1
0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 48 "To see what we got, let \+
Maple plot the function " }{TEXT 305 2 "gg" }{TEXT -1 27 " and the nma
x functions in " }{TEXT 325 2 "cc" }{TEXT -1 79 " (for time t=0 - you \+
can use the subs() command as we did in our example above)" }{TEXT
313 1 "." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 ""
0 "" {TEXT -1 101 "Now we can look how the heat distribution is evolvi
ng. To visualize this, use the animate() command. " }}{PARA 0 "" 0 ""
{TEXT -1 73 "First, create an animation containing gg(x,t) and, for co
mparison, g(x)." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG
{PARA 0 "" 0 "" {TEXT -1 27 "Now add the nmax functions " }{TEXT 312
2 "cc" }{TEXT -1 140 " to your animation as well. Once again, you can \+
see that the high frequency components decay much faster than the low \+
frequency components. " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}
{SECT 1 {PARA 4 "" 0 "" {TEXT -1 8 "Solution" }}{EXCHG {PARA 0 "" 0 "
" {TEXT -1 42 "The Fourier coefficients in a list named " }{TEXT 302
2 "c:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 40 "c := [seq(int(g(x)*f,x=0..
2*Pi),f=fkn)];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 48 "So these are th
e terms in our fourier expansion:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0
42 "cc := [seq(c[i]*fknt[i],i=1..nops(fknt))];" }}}{EXCHG {PARA 0 ""
0 "" {TEXT -1 59 "The approximation of our initial distribution theref
ore is " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "gg := sum(cc[i],i=1..nop
s(fknt));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 14 "And the plots:" }}
{PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "plot([g(x),subs(t=0,gg),seq(subs(t=
0,cc[i]),i=1..nops(fknt))],x=0..2*Pi);" }}}{EXCHG {PARA 0 "" 0 ""
{TEXT -1 48 "Now to the animations. First without components:" }
{MPLTEXT 1 0 0 "" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 68 "with(plots):ani
mate(\{g(x),gg\},x=0..2*Pi,t=0..2,frames=32,color=red);" }}{PARA 0 ""
0 "" {TEXT -1 32 "And with the fourier components:" }}{PARA 0 "> " 0 "
" {MPLTEXT 1 0 75 "with(plots):animate(\{g(x),gg,op(cc)\},x=0..2*Pi,t=
0..2,frames=32,color=red);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0
0 "" }}}}}}{MARK "1 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }
{PAGENUMBERS 0 1 2 33 1 1 }