{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 }