{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 "2D Comment" 2 18 "" 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 } {CSTYLE "" 18 256 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" 18 257 "" 0 1 0 0 255 1 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 258 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 259 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 260 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 } {CSTYLE "" -1 261 "" 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 262 "" 0 1 255 0 0 1 0 0 0 0 0 0 0 0 0 0 }{CSTYLE "" -1 263 "" 0 1 255 0 0 1 0 0 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 "Title" -1 18 1 {CSTYLE "" -1 -1 "Tim es" 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 } } {SECT 0 {EXCHG {PARA 18 "" 0 "" {TEXT -1 42 "Numerical Treatment of Ma thematical Models" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {SECT 1 {PARA 3 "" 0 "" {TEXT -1 25 "Numerical Solution of ODE" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "restart;" }}{PARA 0 "" 0 "" {TEXT -1 51 "We use some packages: we are already familiar with " } {TEXT 0 7 "DEtools" }{TEXT -1 5 " and " }{TEXT 0 5 "plots" }{TEXT -1 1 ":" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "with(DEtools):" }}{PARA 0 " > " 0 "" {MPLTEXT 1 0 12 "with(plots):" }}{PARA 0 "" 0 "" {TEXT -1 104 "The package plottools contains more functions for plotting, espec ially a function to draw a single line:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "with(plottools):" }}{PARA 0 "" 0 "" {TEXT -1 16 "And finally t he " }{TEXT 0 13 "LinearAlgebra" }{TEXT -1 100 "-package - you may gue ss, what kind of functions it contains. In this worksheet, we make use of the " }{TEXT 0 6 "Vector" }{TEXT -1 5 " and " }{TEXT 0 6 "Matrix" }{TEXT -1 176 " data structures defined there. Unfortunately, this pac kage is available only for newer versions of Maple - if you use an old version, you can use lists or the data structures " }{TEXT 0 6 "vecto r" }{TEXT -1 5 " and " }{TEXT 0 6 "matrix" }{TEXT -1 10 " from the " } {TEXT 0 6 "linalg" }{TEXT -1 17 "-package instead." }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 20 "with(LinearAlgebra):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 26 "The populati on model again" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 94 "In the previous \+ worksheet, we solved an ordinary differential equation for a populatio n model:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 50 "popeq := diff(p(t),t) = (0.03 - 0.0005*p(t))*p(t);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "p0 := 2;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "dsolve(\{pop eq, p(0)=p0\},p(t));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "sol := rhs(%);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 60 "DEplot(\{pop eq\},[p(t)],t=0..300,[[p(0)=p0]],linecolor=GREEN);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 166 "As we have a formula for the solution, we can ea sily compute the solution for some value of t - say, TEnd=100. We stor e this value and the corresponding DEplot-graph:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "TEnd := 100;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "pend := evalf(subs(t=TEnd, sol));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "graph_exact := DEplot(\{popeq\},[p(t)],t=0..TEnd,[[p( 0)=p0]],linecolor=GREEN):" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 21 "display(graph_exact);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 4 " " 0 "" {TEXT -1 38 "Preparation for the numerical solution" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 159 "Although we get an exact solution for th is problem, we will solve it again numerically (so we can compare the \+ resulting approximation with the exact solution)." }}{PARA 0 "" 0 "" {TEXT -1 57 "We will compute an approximate solution for the value of \+ " }{XPPEDIT 18 0 "pend = p(TEnd);" "6#/%%pendG-%\"pG6#%%TEndG" }{TEXT -1 12 " from above." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 158 "First, for some integer number n, we div ide the Interval [0,Tend] into n subintervals - here, we choose all su bintervals of the same length delta_t = Tend/n: " }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 7 "n := 5;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 18 "delta_t := TEnd/n;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 46 "The numerical integration will compute values " }{XPPEDIT 18 0 "pnum[k]; " "6#&%%pnumG6#%\"kG" }{TEXT -1 15 ", k=1..n+1 for " }{XPPEDIT 18 0 "p (tp[k]);" "6#-%\"pG6#&%#tpG6#%\"kG" }{TEXT -1 6 " with " }{XPPEDIT 18 0 "tp[k] = (k-1)*delta_t;" "6#/&%#tpG6#%\"kG*&,&F'\"\"\"F*!\"\"F*%(del ta_tGF*" }{TEXT -1 27 ". We store the grid points " }{XPPEDIT 18 0 "tp ;" "6#%#tpG" }{TEXT -1 15 " into a vector:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 36 "tp := ;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 67 "We can draw the intervals with respect to t in the direction field:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 51 "grid \+ := seq(line([tp[k],0],[tp[k],pend]),k=1..n+1):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "display([graph_exact,grid]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 174 "Create \+ an vector with n+1 elements p[1]...p[n+1] (the method itself will not \+ require this storage, but we save all intermediate results, so that we can plot them afterwards):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "p := Vector(1..n+1);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 56 "We know the f irst entry of p from the initial condition:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "p[1] := p0;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 14 "Euler's method" }}{PARA 0 "" 0 "" {TEXT -1 28 "To compute the other values " }{XPPEDIT 18 0 "p[k];" "6#&%\"pG6 #%\"kG" }{TEXT -1 53 ", k=2..n+1, we use Euler's method, that approxim ates " }{XPPEDIT 18 0 "p[k+1];" "6#&%\"pG6#,&%\"kG\"\"\"F(F(" }{TEXT -1 4 " by " }{XPPEDIT 18 0 "p[k]+delta_t*f(t[k],p[k]);" "6#,&&%\"pG6#% \"kG\"\"\"*&%(delta_tGF(-%\"fG6$&%\"tG6#F'&F%6#F'F(F(" }{TEXT -1 8 ", \+ where " }{XPPEDIT 18 0 "f(t,p);" "6#-%\"fG6$%\"tG%\"pG" }{TEXT -1 53 " is the right hand side of the differential equation:" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 88 "for k from 1 to n do\n ft := subs(p(t)=p [k],rhs(popeq));\n p[k+1] := p[k]+delta_t*ft\nod;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 39 "points := [seq([tp[k],p[k]],k=1..n+1)];" }} }{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "pointplot(points,style=LINE ,thickness=2,color=BLUE);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 30 "display([%,graph_exact,grid]);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {SECT 1 {PARA 4 "" 0 "" {TEXT -1 36 "Exercise: Accuracy of Euler's met hod" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 79 "To study the influence of n to the accuracy of the solution, write a procedure " }{TEXT 0 9 "eule rp(n)" }{TEXT -1 36 " that computes the approximation to " }{TEXT 0 4 "pend" }{TEXT -1 53 " for different values of n. Therefore, the result of " }{TEXT 0 9 "eulerp(4)" }{TEXT -1 11 " should be " }{XPPEDIT 256 0 "16.05748735;" "6#-%&FloatG6$\"+N([dg\"!\")" }{TEXT -1 16 ", the res ult of " }{TEXT 0 10 "eulerp(20)" }{TEXT -1 11 " should be " } {XPPEDIT 257 0 "22.44722956;" "6#-%&FloatG6$\"+cHsWA!\")" }{TEXT -1 90 ". You don't need to store the whole array p, it is sufficient to k eep the actual value of " }{XPPEDIT 18 0 "p[k];" "6#&%\"pG6#%\"kG" } {TEXT -1 31 " in a simple (scalar) variable." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "n := 'n';" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 158 "euler p := proc(n)\n local delta_t,i,pnum;\n delta_t := TEnd/n;\n pnum := p0;\n for i to n do \n pnum := pnum + delta_t*subs(p(t)=pnum,rhs( popeq));\n od\nend;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 10 "eul erp(4);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 11 "eulerp(20);" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 5 "pend;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "For the following computations, tell Maple to use 10 digits for floating point computations:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "Digits := 10;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 65 "Then, compute a list of twelve approximations with n=2,4,...2^12:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "res := [seq(eulerp(2^i),i=1..12)]; " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "Plot the result (" }{TEXT 0 11 "pointplot()" }{TEXT -1 35 ", x-axis: 1..12, y-axis: result of " } {TEXT 0 8 "eulerp()" }{TEXT -1 2 "):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "pointplot([seq([i,res[i]],i=1..12)]);" }}}{EXCHG {PARA 0 "" 0 " " {TEXT -1 100 "Next, compute for each of the 12 values the (absolute) error (the difference to the exact solution):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 34 "err := [seq(pend-res[i],i=1..12)];" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 54 "... and plot it in the same way as the previou s plot:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 37 "pointplot([seq([i,err[i] ],i=1..12)]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 111 "To get a better impression of what is happening, plot it again with a logarithmic sca le for the vertical axis (" }{TEXT 0 9 "logplot()" }{TEXT -1 4 ")..." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 47 "logplot([seq([i,err[i]],i=1..12)] ,style=POINT);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 45 ".. and compute \+ and plot the ratio (error for " }{XPPEDIT 18 0 "n = 2^k;" "6#/%\"nG)\" \"#%\"kG" }{TEXT -1 13 ")/(error for " }{XPPEDIT 18 0 "n = 2^(k+1);" " 6#/%\"nG)\"\"#,&%\"kG\"\"\"F)F)" }{TEXT -1 14 ") for k=1..11:" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 46 "pointplot([seq([i,err[i]/err[i+1]], i=1..11)]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 79 "Explain the last r esult with respect to the convergence rate of Euler's method!" }} {PARA 0 "" 0 "" {TEXT 262 307 "Computing err[i+1] - compared with the \+ computation of err[i] - we use twice as many grid points, i.e. divide \+ delta_t by two and receive an error that is (for small values of delta _t) about half the size as err[i] (err[i]/err[i+1] converges towards 2 ): the error is of the same order of magnitude as delta_t!" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 139 "No w, do the same computations (except for the last two plots) again for \+ an floating point accuracy of only 4 digits and explain the result!" } }{PARA 0 "" 0 "" {TEXT 263 182 "Here, the error increases again, if de lta_t becomes too small: the accumulated round-off-error form more and more floating point operations dominates the overall error of the met hod!" }}{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 49 "Numerical solution of the heat t ransport equation" }}{PARA 0 "" 0 "" {TEXT -1 90 "As an example for th e numerical solution of PDEs, we use again a heat transport equation: \+ " }}{PARA 0 "" 0 "" {TEXT -1 6 "solve " }{XPPEDIT 18 0 "diff(u(x,t),t) = diff(u(x,t),`$`(x,2));" "6#/-%%diffG6$-%\"uG6$%\"xG%\"tGF+-F%6$-F(6 $F*F+-%\"$G6$F*\"\"#" }{TEXT -1 23 " for x in the interval " } {XPPEDIT 18 0 "[0, 2*Pi];" "6#7$\"\"!*&\"\"#\"\"\"%#PiGF'" }{TEXT -1 23 " and t in the interval " }{XPPEDIT 18 0 "[0, Pi/2];" "6#7$\"\"!*&% #PiG\"\"\"\"\"#!\"\"" }{TEXT -1 43 " with the initial temperature dist ribution " }{XPPEDIT 18 0 "u(x,0) = sin(x/2)+sin(x);" "6#/-%\"uG6$%\"x G\"\"!,&-%$sinG6#*&F'\"\"\"\"\"#!\"\"F.-F+6#F'F." }{TEXT -1 29 " and t he boundary conditions " }{XPPEDIT 18 0 "u(0,t) = 0;" "6#/-%\"uG6$\"\" !%\"tGF'" }{TEXT -1 5 " and " }{XPPEDIT 18 0 "u(2*Pi,0) = 0;" "6#/-%\" uG6$*&\"\"#\"\"\"%#PiGF)\"\"!F+" }{TEXT -1 9 " for all " }{XPPEDIT 18 0 "t;" "6#%\"tG" }{TEXT -1 1 "." }}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 38 "Preparation (grids and initial values)" }}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 8 "restart;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "with( LinearAlgebra):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 12 "with(plots):" }} }{EXCHG {PARA 0 "" 0 "" {TEXT -1 22 "Divide the x-interval " } {XPPEDIT 18 0 "[0, 2*Pi];" "6#7$\"\"!*&\"\"#\"\"\"%#PiGF'" }{TEXT -1 6 " into " }{TEXT 0 2 "nx" }{TEXT -1 30 " subintervals of equal length " }{TEXT 0 2 "hx" }{TEXT -1 41 " and store the grid points in the vec tor " }{TEXT 0 2 "xp" }{TEXT -1 1 ":" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 8 "nx := 5;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "hx := 2*Pi/nx;" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "xp := Vector(nx-1,i->i*hx);" }}} {EXCHG {PARA 0 "" 0 "" {TEXT -1 31 "Do the same for the t-interval " } {XPPEDIT 18 0 "[0, Pi/2];" "6#7$\"\"!*&%#PiG\"\"\"\"\"#!\"\"" }{TEXT -1 2 " (" }{TEXT 0 5 "nt=40" }{TEXT -1 15 " subintervals):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 13 "TEnd := Pi/2;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 9 "nt := 40;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 14 "ht := \+ TEnd/nt;" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "tp := Vector(nt+1,j->(j -1)*ht);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 37 "The initial temperatu re distribution:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 25 "u0 := x->sin(x/ 2)+sin(x);" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 44 "u0plot := plot(u0(x), x=0..2*Pi,thickness=2):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "display( u0plot);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 57 "For the numerical sol ution, we store for each grid point " }{TEXT 0 5 "xp[i]" }{TEXT -1 14 " the value of " }{TEXT 0 9 "u0(xp[i])" }{TEXT -1 1 ":" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 39 "u0num := < seq(u0(xp[i]), i=1..nx-1) >;" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 82 "u0numplot := pointplot([seq( [xp[i],u0num[i]], i=1..nx-1)],style=LINE,thickness=2):" }}{PARA 0 "> \+ " 0 "" {MPLTEXT 1 0 28 "display([u0numplot,u0plot]);" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 36 "Finite difference approximation for " }{XPPEDIT 18 0 "diff(u(x, t),`$`(x,2));" "6#-%%diffG6$-%\"uG6$%\"xG%\"tG-%\"$G6$F)\"\"#" }{TEXT -1 0 "" }}{PARA 0 "" 0 "" {TEXT -1 41 "As a finite difference approxim ation for " }{XPPEDIT 18 0 "diff(u(x,t),`$`(x,2));" "6#-%%diffG6$-%\"u G6$%\"xG%\"tG-%\"$G6$F)\"\"#" }{TEXT -1 9 ", we use " }{XPPEDIT 18 0 " (u(x+hx,t)-2*u(x,t)+u(x-hx,t))/(hx^2);" "6#*&,(-%\"uG6$,&%\"xG\"\"\"%# hxGF*%\"tGF**&\"\"#F*-F&6$F)F,F*!\"\"-F&6$,&F)F*F+F1F,F*F**$F+F.F1" } {TEXT -1 71 ". We store the linear operator that maps a vector contain ing u-values (" }{TEXT 0 5 "u0num" }{TEXT -1 163 ", e.g.) to the appro ximation of the second derivative in a matrix A. From the viewpoint of efficiency, this isn't a good idea, but it keeps our maple code small ..." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "A := BandMatrix([[(1/ hx^2)$(nx-2)],[(-2/hx^2)$(nx-1)],[(1/hx^2)$(nx-2)]]);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 125 "You may ask about the first and the last grid point - they are handled correct by the above discretization as \+ in our example " }{XPPEDIT 18 0 "u(0,t);" "6#-%\"uG6$\"\"!%\"tG" } {TEXT -1 5 " and " }{XPPEDIT 18 0 "u(2*Pi,t);" "6#-%\"uG6$*&\"\"#\"\" \"%#PiGF(%\"tG" }{TEXT -1 14 " are always 0." }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 28 "Now, we can ' differentiate' " }{TEXT 0 5 "u0num" }{TEXT -1 1 ":" }}{PARA 0 "> " 0 " " {MPLTEXT 1 0 16 "\014Au := A. u0num;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 27 "The result - compared with " }{XPPEDIT 18 0 "diff(u0(x),` $`(x,2));" "6#-%%diffG6$-%#u0G6#%\"xG-%\"$G6$F)\"\"#" }{TEXT -1 20 " - looks reasonable:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 76 "Auplot := poi ntplot([seq([xp[i],Au[i]], i=1..nx-1)],style=LINE,thickness=2):" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 55 "dduplot := plot(diff(u0(x),x,x),x=0 ..2*Pi,thickness=2):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "display([Au plot,dduplot]); " }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {SECT 1 {PARA 3 "" 0 "" {TEXT -1 38 "Time stepping: (explicit) Euler s cheme" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 133 "Now, we can perform an e xplicit Euler method with respect to t as in the case of the ODE above : we approximate the left hand side of " }{XPPEDIT 18 0 "diff(u(x,t),t ) = diff(u(x,t),`$`(x,2));" "6#/-%%diffG6$-%\"uG6$%\"xG%\"tGF+-F%6$-F( 6$F*F+-%\"$G6$F*\"\"#" }{TEXT -1 4 " by " }{XPPEDIT 18 0 "(unum[x,t+ht ]-unum[x,t])/ht;" "6#*&,&&%%unumG6$%\"xG,&%\"tG\"\"\"%#htGF+F+&F&6$F(F *!\"\"F+F,F/" }{TEXT -1 28 " and the right hand side by " }{XPPEDIT 18 0 "A*unum;" "6#*&%\"AG\"\"\"%%unumGF%" }{TEXT -1 2 ". " }}{PARA 0 " " 0 "" {TEXT -1 47 "The function timestep() computes the values of " } {TEXT 259 4 "unum" }{TEXT -1 38 " in the next time step t+ht for given " }{TEXT 258 4 "unum" }{TEXT -1 11 " at time t:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 35 "timestep := u -> evalf(u + A.u*ht);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 18 "If we apply it to " }{TEXT 0 5 "u0num" }{TEXT -1 38 ", we get an approximation for u(x,ht):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 16 "timestep(u0num);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 129 "To plot the solution, we store all values in a matrix u3d - each \+ column corresponds with one grid point with respect to the time " } {TEXT 260 1 "t" }{TEXT -1 63 ", while the rows correspond to the grid \+ points with respect to " }{TEXT 261 2 "x." }}{PARA 0 "" 0 "" {TEXT -1 55 "Therefore, the first column contains the initial values" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "u3d := evalf(u0num);" }}}{EXCHG {PARA 0 " " 0 "" {TEXT -1 18 "Do the first step:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "unum := timestep(u0num);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 37 "Add it as a new column to the matrix:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "u3d := ;" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 23 "And do all other steps:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 74 "for i from 2 to nt do \n unum := timestep(unum);\n u3d := \nod:" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 19 "This is the resu lt:" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 27 "matrixplot(u3d,axes=BOXED); " }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 17 "Try it again for " }{TEXT 0 5 "nx=20" }{TEXT -1 5 " and " }{TEXT 0 5 "nx=50" }{TEXT -1 1 "!" }} {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 44 "Exercise: co mparison with the exact solution" }}{EXCHG {PARA 0 "" 0 "" {TEXT -1 92 "Your job is to compare the numerically computed values to the exac t solution of the problem." }}{PARA 0 "" 0 "" {TEXT -1 25 "First, defi ne a function " }{TEXT 0 8 "uex(x,t)" }{TEXT -1 108 " that computes th e solution of the initial value problem (we know the solution from the previous worksheet):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 52 "uex := (x, t) -> sin(x/2)*exp(-t/4) + sin(x)*exp(-t);" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 22 "Next, create a matrix " }{TEXT 0 5 "uex3d" }{TEXT -1 23 " with the same size as " }{TEXT 0 3 "u3d" }{TEXT -1 118 " and store th e values of uex at the grid points used in the Euler method. Store the m in the same order as we used for " }{TEXT 0 3 "u3d" }{TEXT -1 1 "!" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 58 "uex3d := Matrix(nx-1,nt+1,(i,j)-> evalf(uex(xp[i],tp[j])));" }}}{EXCHG {PARA 0 "" 0 "" {TEXT -1 28 "Plot the difference between " }{TEXT 0 5 "uex3d" }{TEXT -1 5 " and " } {TEXT 0 3 "u3d" }{TEXT -1 6 " for " }{TEXT 0 4 "nx=5" }{TEXT -1 4 "an d " }{TEXT 0 5 "nx=20" }{TEXT -1 1 "!" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 33 "matrixplot(uex3d-u3d,axes=BOXED);" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 3 "" 0 "" {TEXT -1 58 "Additional e xercise (not mandatory): implicit Euler scheme" }}{PARA 0 "" 0 "" {TEXT -1 193 "If you want to get more experience with the solution of \+ PDE of this type, you may want to see as scheme that is unconditionall y stable (i.e. for all values of hx and ht without 'exploding' for " } {XPPEDIT 18 0 "hx^2/2 < ht;" "6#2*&%#hxG\"\"#F&!\"\"%#htG" }{TEXT -1 99 "). The implicit Euler scheme has this desired feature, but for the cost of more computational work." }}{PARA 0 "" 0 "" {TEXT -1 50 "The \+ key idea is to discretize the right hand side " }{XPPEDIT 18 0 "diff(u (x,t),`$`(x,2));" "6#-%%diffG6$-%\"uG6$%\"xG%\"tG-%\"$G6$F)\"\"#" } {TEXT -1 15 " of the PDE by " }{XPPEDIT 18 0 "A*unum[x,t+ht];" "6#*&% \"AG\"\"\"&%%unumG6$%\"xG,&%\"tGF%%#htGF%F%" }{TEXT -1 13 " (instead b y " }{XPPEDIT 18 0 "A*unum[x,t];" "6#*&%\"AG\"\"\"&%%unumG6$%\"xG%\"tG F%" }{TEXT -1 33 " from the explicit Euler scheme)." }}{PARA 0 "" 0 " " {TEXT -1 42 "Then, we get a system of linear equations " }{XPPEDIT 18 0 "(I-ht*A)*unum[t+ht] = unum[t];" "6#/*&,&%\"IG\"\"\"*&%#htGF'%\"A GF'!\"\"F'&%%unumG6#,&%\"tGF'F)F'F'&F-6#F0" }{TEXT -1 8 ", where " } {XPPEDIT 18 0 "I;" "6#%\"IG" }{TEXT -1 49 " denotes the identity matri x of the same size as " }{XPPEDIT 18 0 "A;" "6#%\"AG" }{TEXT -1 1 "." }}{PARA 0 "" 0 "" {TEXT -1 23 "Use the Maple function " }{TEXT 0 13 "L inearSolve()" }{TEXT -1 39 " to write a new time stepping function " } {TEXT 0 14 "timestep_imp()" }{TEXT -1 29 " that implements this method :" }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 70 "timestep_imp := u -> ev alf(LinearSolve(IdentityMatrix(nx-1)-ht*A, u));" }}}{PARA 0 "" 0 "" {TEXT -1 229 "To test your function and to compare it to timestep() wi thout waiting to long for the results, you should use an example with \+ small nx=5 and nt=40 and produce instability of the explicit method by choosing a large value for TEnd=" }{XPPEDIT 18 0 "50*Pi/2;" "6#*(\"#] \"\"\"%#PiGF%\"\"#!\"\"" }{TEXT -1 1 "." }}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 20 "timestep_imp(u0num);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "u3d_imp := evalf(u0num):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "unum := timestep_imp(u0num):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "u3d_imp := :" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 86 "for i from 2 to nt do\n unum := timestep_imp(unum); \n u3d_imp := \nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 31 "matrixplot(u3d_imp,axes=BOXED);" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{SECT 1 {PARA 4 "" 0 "" {TEXT -1 24 "Maki ng it more efficient" }}{PARA 0 "" 0 "" {TEXT -1 290 "You may be not s atisfied to solve the linear system again for every time step. We can \+ make things better by storing the Cholesky decomposition L of the coef ficient matrix (LUDecomposition() does this for us) and tell LinearSol ve to use it (and, therefore, only to do the back substitution):" }} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 73 "L := evalf(LUDecomposition(I dentityMatrix(nx-1)-ht*A,method='Cholesky'));" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 49 "timestep_imp2 := u -> evalf(LinearSolve([L], u)) ;" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 24 "u3d_imp := evalf(u0num ):" }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 29 "unum := timestep_imp2(u0num): " }}{PARA 0 "> " 0 "" {MPLTEXT 1 0 28 "u3d_imp := :" } }{PARA 0 "> " 0 "" {MPLTEXT 1 0 87 "for i from 2 to nt do\n unum := t imestep_imp2(unum);\n u3d_imp := \nend:" }}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}}{EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}} {EXCHG {PARA 0 "> " 0 "" {MPLTEXT 1 0 0 "" }}}{EXCHG {PARA 0 "> " 0 " " {MPLTEXT 1 0 0 "" }}}}{MARK "3 0 0" 0 }{VIEWOPTS 1 1 0 1 1 1803 1 1 1 1 }{PAGENUMBERS 0 1 2 33 1 1 }