{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 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 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 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 47 "Partial differential equ ations: heat conduction" }}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 36 "The on e dimensional heat conduction " }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 295 "We want to consider the problem of heat conducting in a medium (witho ut currents or radiation) in the one dimensional case. We are given a \+ wire (of infinite length) which has a given distribution of temperatur e at time t=0. The question is how the heat is conducted through the b ody of the 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 "heat equation" }{TEXT -1 144 ". It is an example of a partial \+ differential 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*dif f(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. The heat equation shows what a PDE is in general: We are g iven some kind of relationship between a function u(x,t) and/or its pa rtial derivatives. In the one dimensional case this can be partial der ivatives in time t or in space x. " }}{PARA 0 "" 0 "" {TEXT -1 172 "Th e problem is to find functions u(x,t) which satisfy the given relation ship and sometimes satisfy additional conditions (initital value condi tions or boundary conditions)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }} {PARA 0 "" 0 "" {TEXT -1 270 "PDEs like the heat equation (but usually much more complex) very often arise in modeling problems in science, \+ engineering and financial engineering. Sometimes we are able to find e xact solutions using analytic methods. Maple is also capable of solvin g some PDEs via the " }{TEXT 277 11 "pdesolve( )" }{TEXT -1 16 " comma nd (Click " }{HYPERLNK 17 "here" 2 "pdesolve" "" }{TEXT -1 199 " for d etails of this command). However, often a PDE cannot be solved analyti cally and our only chance is a numerical (and therefore approximate) s olution (we will discuss that later in this lecture)." }}{PARA 0 "" 0 "" {TEXT -1 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 132 "In this exerc ise we want to solve the heat equation exactly. Older Versions of Mapl e cannot solve 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 techniq ue to solve this type of equations without Maple is " }{TEXT 278 13 "b y 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 "As X(x) is independent of t and T(t) of x the partial de rivatives become much easier." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "se p := eq/X(x)/T(t);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 290 "We now hav e \"separated\" the variables: each side of the equation depends only \+ on one variable. The only way this equation can be satisfied is if bot h sides are constant. We only consider the case of negative values (th is is the case that is of physical interest) and denote the constant b y " }{XPPEDIT 18 0 "-omega^2;" "6#,$*$%&omegaG\"\"#!\"\"" }{TEXT -1 53 ". So 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 \+ equation (as T(t) is a function of one variable). This equation can ea sily 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 "The same can be done for the right hand side of the equat ion:" }}{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 h ave the constant (due to integration) " }{TEXT 283 3 "_C1" }{TEXT -1 8 " and in " }{XPPEDIT 18 0 "Xsol;" "6#%%XsolG" }{TEXT -1 15 " the con stants " }{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 "C 1 " }{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 con stant due to our separation ansatz):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 32 "Xsol := 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 "om ega;" "6#%&omegaG" }{TEXT 288 15 ", _C1, _C2, _C3" }{TEXT -1 339 " can be chosen arbitrarily. Usually we are given some initial value condit ions e.g. the values u(x,t=0) which would in our case be the temperatu re distribution at t=0. Using this additional information we can calcu late the constants and obtain a unique solution. In this exercise we j ust 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,usol);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 29 "Now we can plot th e solution." }}{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 resul t to 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 "" }}}}}} {MARK "0 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }