Introduction to Applied Mathematics
Introduction to Applied Mathematics MA 325
Popular in Course
Chem 0014A (Scerri)
verified elite notetaker
Popular in Mathematics (M)
verified elite notetaker
This 18 page Class Notes was uploaded by Braeden Lind on Thursday October 15, 2015. The Class Notes belongs to MA 325 at North Carolina State University taught by Robert White in Fall. Since its upload, it has received 29 views. For similar materials see /class/223679/ma-325-north-carolina-state-university in Mathematics (M) at North Carolina State University.
Reviews for Introduction to Applied Mathematics
Report this Material
What is Karma?
Karma is the currency of StudySoup.
Date Created: 10/15/15
Lecture 2 Discrete Versus Continuous Models Introduction Initial value problems have the form ut ftu and u0 given 1 The simplest cases can be solved by separation of variable but in general the closed form solutions are dif cult to nd for example ftu t2 uz Therefore one is forced to consider various approximation methods In this lecture we consider the Euler nite difference method This method is based upon approximating the initial value problems by a sequence of calculations which divide the time interval 0T into K parts The mesh size is h TK We shall see that the error made by the Euler approximation is proportional to h Again consider a wellstirred liquid such as a cup of coffee We will assume that the temperature is uniform with respect to space but the temperature may be changing with respect to time We wish to predict the temperature as a function of time given some initial observations about the temperature For example a cup of coffee will cool much faster in a metal cup than in an insulated cup Or the coffee will cool much faster if the surrounding temperature is very cold Model The continuous version of Newton39s empirical law of cooling states that the rate of change of the temperature is proportional to the difference in the surrounding temperature and the temperature of the liquid ut cusur u where u is the temperature of the liquid u0 is observed usur is the surrounding temperature and c is the coefficient and re ects the insulation of the container If c 0 then there is perfect insulation and the liquid temperature must remain at its initial value For large c the liquids temperature will rapidly approach the surrounding temperature The closed form solution of this differential equation can be found by the separation of variables method and is for usur equal a constant ut usur 110 11239 t If the c is not given then it can be approximated by using a second observation ut1 111 In this case one can approximate c by a nite difference approximation which we will later illustrate If usur is a function of t a closed form can be found solution provided it is not too complicated to be easily integrated Method Euler s method involves the approximation of ut by the finite difference uk1 ukh where h TK uk is an approximation of ukh and f is evaluated at khuk If T is not nite then h will be fixed and k may range over all of the positive integers The differential equation 1 can be replaced by uk ukh fk1huk1 or uk1ukh fkhuk fk1huk12 or uk ukh fkhuk 2 The third choice is the simplest because it does not require the solution of a possibly nonlinear problem at each time step The scheme given by 2 is call Euler39s method and it is a discrete model of 1 Example Let ftu t2 u2 and u0 0 Suppose we let h l and wish to approximate ul Then we must compute u10 by repeating the calculation in the loop 10 times So u0 uO 0 u1 1 102 12 11 m uattime 1 and u2 11112 112 1222 m u attime 2 What is the appropriate choice for h Implementation Use Newton39s law of cooling to predict the temperature of a cup of coffee The coffee is kept well stirred by walking around the classroom Also we need to know when its temperature reaches 110 after that the professor gets irritable Today the initial temperate was 200 and five minutes later it was observed to be 190 The second measurement varies with the type of coffee container The continuum model is ut cusur u where u0 200 u5 190 and usur 70 It is possible to solve this model in closed form and the solution is ut 70 200 70equot The c is determined from the second observation and is c l5ln12l3 01600854 The discrete model is formed by using Euler39s algorithm and is uk1 uk hcusur uk where u0 200 um 70 T 50 h 50K The c must be determined from the second observation Here the increase of time from 0 to 5 is not too large and so we may approximate the ut0 by 190 2005 and now put this into the differential equation 190 2005 c70 200 or c70 190 or c70 200 1902 We will choose the last possible right side so that c 2125 016 which is not too far from the c in the continuum model For K 10 the time step is h 5010 The computations below use a larger c 061 and indicate the effect of the mesh size on the numerical solution For this example the numerical solution appears to converge as K increases from 4 to 8 and to 16 The graphs of the numerical solution are in increasing order The next numerical experiments are done using Matlab They illustrate numerical error as a function of the number of steps This done with the equation for Newton39s law of cooling with usur 70 and c 061 Matlab has a number of solvers for ordinary differential equations Matlab Code eulerdem y1 200 T KK 16 h TKK t1 0 for k 1KK yk1 yltkgt hfeu1tkyk tk1 tk h end p10tty Matlab Function feulm function feul feultx feul 06170 x EDU eulerde KK 4 EDU t4 t EDU y4 39 EDU eulerde K 8 EDU t8 t39 EDU y8 y EDU eulerde KK 16 EDU t16 t39 EDU y16 y EDU whitebg EDU plott4y4t8y8t16y16 Figure Convergence of Numerical Solution In the following table We list the discretization errors for the above rst order Euler algorithm and the second order EulerTrapezoid algorithm The discretjzatjon errors are the difference in the solution to differential equation and the numerical method with no roundoff errors The tabulated errors Were the largest values of all the points in time The constant c Was from the discrete model Table Errors K Euler Error Euler Ermth ETmp Err ETrap Err112 10 19710 3942 025600 001 20 09664 3866 006400 001 40 04786 3829 001600 001 80 02382 3811 000399 001 160 01188 3802 000010 001 Assessment If the liquid is not well stirred then the temperature will vary within the container This results in diffusion of heat from a warm region to a cool region The above model does not take this into consideration Also in our analysis the surrounding temperature was held constant Variable surrounding temperature can easily be accounted for in our present model Clearly this is a very inexpensive algorithm to execute However there are sizable errors There are two types of errors Discretization error E uk ukh where uk is from Euler39s algorithm 2 with no roundoff error and ukh is from the exact continuum solution 1 Roundoff error Ef Uk uk where Uk is from Euler39s algorithm but with oating point numbers The overall error contains both errors and is E Ef Uk ukh In the above table the discretization error for Euler s method is bounded by a constant times h and the discretization error for the EulerTrapezoid method see the line above 2 is bounded by an constant times h squared Homework 1 Write a computer code for Euler39s algorithm use the above example with ftu t2 uz Experiment with different initial conditions uO and different step lengths h 2 Verify the above calculations for the discretization error 3 Consider the above coffee cup problem but now the classroom temperature is observed to increase in a linear way from 70 to 80 at the end of class 50 minutes later Find both the continuum and discrete solutions and do the computations similar to those given in the above table Lecture 5 Analysis of Stability Introduction In this lecture we continue to consider heat conduction in a thin electrical wire which is thermally insulated on its surface The time dependent model of the temperature will have the form newy A oldy b In general the matrix A can be extremely large but it will also have a special structure with many more zeros than nonzero components As time marches one can expect the temperature to reach a steady state so that y A y b In the previous lecture we formulated the discrete model based on the Fourier law for heat diffusion in a thin wire Explicit Finite Difference Model of Heat Transfer uik1 Atpc f 0cui1k ui1k l 20c uik where l i lnl k 0maxkl uio 0 fori lnl and 2 u0kunk 0 fork lmaxk 3 Equation 2 is the initial temperature set equal to zero and 3 is the temperature at the left and right ends set equal to zero Equation 1 may be put into the matrix version of the first order finite difference method For example if the wire is divided into four equal parts then n 4 and 1 may be written as either as three scalar equations or as one 3D vector equation u1k1 Atpc f 0qu 0 l 20cu1k uzk1 Atpc f 0cu3k ulk l 20c uzk 113 1 Atpc f 0c0 uzk 1 20c 113k uk1 A uk b where uf 1 1 20 at uk u bAtpcfl andA a 1 20 a 4 u l a 1 20 An extremely important restriction on the time step At is required to make sure the algorithm is stable Here a l 20c gt 0 is needed to avoid a blowup of the numerical solution This simple condition implies that the matrix products Ak will converge to the zero matrix This will make sure there are no blowups provided the source term f is bounded Stability Condition for 1 l20cgt0and0cgt0 Matrix Products In order to carefully study stability we will need some basic facts about matrices The simplest matrices also called arrays are row and column vectors For example A l 3 5 is a 3x1 row vector and 8 B 1 is a 4x1 column vector 10 If the row and column vectors have the same number of components then we de ne their product to be a single number given by the sum of the products of the components For example A1 4 0 8 2 B 9 10 AB124709810110 De nition Let A ak T denote a lxn row vector where k varies from 1 to n Let B bk denote a nxl column vector where k varies from 1 to n De ne the row column matrix product of A and B to be a single number AB akbk k1 If a matrix has m rows and n columns then we may denote it by A aik where ivaries from 1 to m andj varies from 1 to n In this case A is called an mxn matrix IfA is mxn and B bki is nxp so thatj varies from 1 to p then we can define a matrix product AB to be a mxp matrix whose ij components are row i ofA times columnj of B For example l 165 89 234 78 A B 6 10 11 12 1 5 9 156559 88 100 112 e 718599 152 176 200 128 152 176 20039 76 88 100 112 It is important to note that the product BA may not be de ned even if AB is defined In general BA and AB are not equal for matrices A and B De nition Let A aik be mxn and B bki be nxp The matrix product AB is de ned to be a mxp matrix whose ij component is the product of row I of A and column j of B that is AB 2611ka k1 De nition Let B and D both be nxp matrices The matrix sum of B and D is de ned to be a third nxp matrix whose kj components are the sum of the kj components of B and D that is BD bk dkj Proposition Let A be mxn B and D be nxp and C be pxq matrices Then the following are true 1 ABD AB AD distributive property 2 ABC ABC associative property Proof of 1 BD bk dkj AB D 16114ka dkj i 611ka alk dkj by the distributive property for real numbers k1 i 611ka 2m alk dkj by the de nition of the summation operation k1 k1 Z 611ka Z 611de by the de nition of matrix addition 161 k1 AB AD by the de nition of matrix products Proof of 2 MC i at Z bye1 F Z alk bkj c by de nition of the summation symbol k1 1 m p Z Z 611ka c by the associative property for real numbers k1 1 p m Z Z all kbkj c 11 by commutative property of addition 1 k1 ABC by the de nition of matrix products Analysis Of Stability The heat conduction in a thin wire has a number of approximations Different mesh sizes in either the time or space variable will give different numerical results However if the stability conditions holds and the mesh sizes decreases then the numerical computations will differ by smaller amounts In the scalar version of rst order nite difference models the scheme was stable when a l 20c gt 0 In this case yk1 converged to the steady state solution y a y b This is also true of the matrix version of 1 provided the stability condition is satisfied In this case the real number a will be replaced by the matrix A and Ak will converge to the zero matrix Steady State Theorem Consider the matrix version of the first order nite difference equation yk1 A yk b where A is a square matrix If Ak converges to the zero matrix and y Ay b then regardless of the initial choice for yo yk converges to y 1 Proof Subtract yk A yk b and y A yk b and use the properties of matrix products to get y 1 7yAykbAyb Ayk y by the distributive property AAyk1 y by recursion A2 yk39l y by the associative property Ak y0 y Since Ak converges to the zero matrix the column vectors yk1 y must converge to zero The assumption that Ak must converge to the zero matrix seems difficult to verify However by using matrix norms we can formulate an easy test which ensures the stability conditions on the time step will imply this assumption is true Homework 1 Consider the 3x3 A matrix in 4 Observe Ak for different values of or so that the stability condition either does or does not hold 2 Vary 4 so that there are n 5 equal parts of the wire 3 Find 2x2 matrices A and B such that AB is not equal to BA Lecture 6 Discrete Fourier Transform In the previous lecture we introduced the discrete Fourier transform as given either by summations or as a matrix vector product The discrete Fourier transform off fg 1 is another vector whose k1h component is r171 r171 gammyn Iannle z j j j0 Or in terms of a matrix product it is Ff where F has kj component 21 Three very important properties include trig functions inverse discrete Fourier transforms and the convolution identity These will allow us to formulate a lter to purge images of unwanted periodic noise see the Matlab code fftsinem The Matlab code ffttrigm computed the discrete Fourier transform for the sine and cosine functions with different frequencies f 100 and f 40 Note the spikes occur at frequencies f and n f The amplitudes are proportional the inZ will be positive real for the cosine function and opposite sign and imaginary for the sine function 4nm l l l 3000 39 realfrtscos2pi100 2000 7 7 10007 e 0 gt f T 1000 l l l l l 200 400 600 800 1000 1200 frequency qnnn l l l l l 1000 n imagf 4sin2pi40t J V 71000 7 i 7 7000 l l l l l 200 400 600 800 1000 1200 frequency Discrete Fourier Transform of Trig Functions Let t 0 lnll The discrete Fourier transform of en ft is 0 for k i f and n for k f The discrete Fourier transform of e39izmis 0 for k i nf and n for k nf The discrete Fourier transform of cos271f t is 0fork fornf n2 for k fand nf The discrete Fourier transform of sin271f t is 0fork fornf n2i for k fand n2 i for k nf The proofs of the Fourier transforms for the sine and cosine vectors follow from the Fourier transforms of the exponential functions and Euler s formula The above example also illustrates these properties In order to understand the first formula consider the special case n 3 and use the matrix representation for the discrete Fourier transform 1 1 1 em 1 1 1 1 FEMrt 1 z1 z2 emf 3 1 z1 z2 z391f 1 z2 z1 emf 3 1 z2 z1 z392f Iff0 391 1 1quot1 111 3 Fe z m 1 z1 z2 l lzz2 0 1 z2 zll lzzz 0 IffL l l l l lZZZ 0 Few 1 21 22 2 1 111 3 1 z2 zlz392 lzlz2 0 Iff2 391 1 1quot 1 1zz2 0 Fe z m 1 z1 z2 z392 lzzz 0 1 z2 zlz391 lll 3 The reader should verify the remaining Fourier transforms The Inverse Discrete Fourier Transform The inverse discrete Fourier transform of g g0 g1 gn1 is another vector h whose k component 1s 1 quot71 1 quot71 127rkjn l27rVlky Ze g quotXe g 71 F0 71 F0 quot71 Z 0 7k g IH Or in terms of a matrix product it is l F39lg where F1 has 939 component z39k n In order to con rm this again consider the special case n 3 111 lzzz 1zz2 1zz2 111 lzzz lzzz 1zz2 111 3 0 0 1 0 3 0 1 3 0 0 3 The inverse Fourier transform can be efficiently computed by using the same technique as in the fast Fourier transform The following illustrates the Matlab command ifft for computing the inverse Fourier transform gtgt f 1 7 2 gtgt g f g 100000 35000 43301i 35000 43301i gtgt ifftg ansl0000 70000 20000 The Discrete Convolution and Fourier Transform The discrete convolution of two vectors of length n is vector of length 2n1 The components correspond to the 2n2 degree polynomial that is a product of the two nl degree polynomials with coefficients from the two vectors of length n For example for n 3 and the two vectors a a0 a1 a2 and b b1 b1 b2 pa xpb x a0 611x1 612x2 b0 blx1 bzxz czoboaob1 b0alx1 aob2 alb1 azbo x2 61le azblx3 azb2 x The convolution of a and b E convab aob0 aob1 b0al aob2 alb1 azb0 61le azb1 azbz Then me p xpbx The analogue of the convolution property for the discrete Fourier transforms requires the use of padded discrete Fourier transforms This is a technical difficulty caused be the fact that the dimension of vectors increase during the convolution of two vectors Let a be a n vector A padded a vector is a 2n 1 vector with a in the first n components and zeros in the last n 1 components For example let n 3 and a a0 a1 612 A a2 Its Fourier transform is 0 0 l l l l 1 a0 1 l l 1 z z2 z3 24 al 1 z 22 FA 1 z2 z4 z 23 a2 1 610 22 611 24 a2 1 23 z z4 z2 0 1 23 z 1 z4 z3 z2 z 0 1 z4 23 The component wise product with a Fourier transform of a padded vectors of a and b is l l l l l l 1 z z2 1 z z2 FAFB1 610 22 611 z4 a2l b0 z2 b1 24 b2 1 23 z 1 23 z 1 z4 z3 1 z4 z3 l l l 1 z z2 1 a0b0 z2 a0b1a1b0 z4 a0b2a1b1a2b0 l 23 z 1 z4 23 l l 23 z4 3 z a1b2a2b1 z azbz z4 z2 z2 z F conva b Convolution Propelty Let A and B be the padded vectors of two 71 vectors a and b The 271 dimensional Fourier transform of the convolution of a and b is the point wise or array product of the Fourier transforms of A and B Fconvw b FA PB Example Let n 3 and a 2 3 6 and b 8 6 3 gtgt a 2 3 6 gtgt b 8 6 3 gtgt convab ans 16 36 72 45 18 gtgtAa00 gtgtBb00 gtgt ta fftA ffta 110000 63799i 19271 63799i gtgt fftb fftB fftb 170000 74697i 74271 74697i gtgt lhs fftconvab lhs 10e002 18700 03299i 06197 03299i gtgt rhs fftafftb rhs 10e002 18700 03299i 06197 03299i 14271 39430i 14271 39430i 19271 40729 06735i 40729 06735i 74271 00847 01510i 00847 01510i 06197 0084701510i 0084701510i 06197 The convolution property can be restated as conva b 117419 14 FB If B is an image then its Fourier transform will re ect the frequencies of the periodic parts of the image By masking or ltering out the unwanted frequencies one can obtained a new image by applying the inverse Fourier transformation A lter is a matrix with the same dimension as the Fourier transform of the padded image whose components vary from 0 to 1 If the component is 1 then the frequency is allowed to pass if the component is 0 then the frequency is tossed out Let Filter represent such a matrix Then the ltered image is NewB F391Filter FB Application to Low Pass Filter 20 i i i i i i i i i Periodic Noise on Sine Curve i i 0 01 02 03 04 05 06 07 08 09 1 6000 i i i i Shi ed Fourier Transform 4000 2000 e r 0 500 1000 1500 2000 2500 Matlab Code fftsinem t0 Filter High Frequencies Clear t 7 0001l Define 10 sin2 pi t with high frequency variations x 10sin2pit 2sin2pi200t lcos2pi210t fftx fftx2002 Padded fft fftx fftshiftfftx figurel subplot2ll plottx subplot2l2 plotabsfftx freqfilter zerosl2002 Low pass filter freqfilterl700l300 l newfftx freqfilterfftx newx ifftshiftnewfftx newx ifftnewx2002 Padded inverse fft figure2 subplot2ll plotabsnewfftx subplot2l2 plottrealnewxl121001 6000 i 4000 Filtered Frequency Domain 2000 e e i i i 0 500 1000 1500 2000 2500 10 i i Filtered Time Domain 0 01 02 03