Differential Equations in Matlab Cheng Ly1 1

University of Pittsburgh, Department of Mathematics, Pittsburgh, Pennsylvania

15260, USA. E-mail: [email protected]

This workshop assumes you have some familiarity with ordinary (ODEs) and partial differential equations (PDEs), and also that you are relatively comfortable with basic programming in Matlab. Matlab is a very powerful high-level computing tool with a lot of nice built in packages. It is relatively easy to learn, but lags in computation time compared to complied languages such as Fortran, C, or C++. In particular, when solving PDEs with 3 or more spatial dimensions (i.e., 4 or more when including time) Matlab is not something you should use. However, for 2 or less spatial dimensions, Matlab is fast enough these days to be successfully utilized. Numerical simulations are undoubtedly an important tool in investigating the complex behavior of many systems (presumably that is why you are here – unless you only came for the free pizza in which case you should...). Mathematical analysis and ingenuity are obviously important as well, but the value of numerical tools cannot be underscored.

1. Getting Started: A Quick Review Matlab is optimized to work with vectors and matrices. You can specify variables to be row or column vectors. For example, the command >> x=[1 2 3]; creates a 3x1 row vector with the entries 1, 2, and 3. The semicolon ; at the end hides the output (try it without a semicolon). You can also create a column vector with the command >> x=[3; 2; 1]; The semicolons between [ and ] terminates that the row. You can access parts of the vector x. To see how this works, type >> x(1:2) and 1

>> x(2:end) and >> x(1:end-1) You can also reverse the order of the entries with a command like >> x(end:-1:2) Important functions. The function zeros(n,m) and ones(n,m) creates matrices of size n x m with all 0’s and 1 entries (respectively). Also, eye(N) creates the identity matrix of size N x N. You can concatenate vectors and columns easily in between brackets. Type in the following commands: >> A=[1 0 -1; 2*ones(3,1); 8 zeros(2,1)] >> B=[eye(3); A] You can also multiply and divide component-wise when you have a matrix/vector with the . operator. Type: >> x./3 to divide every entry in x by 3. The function diag(x) creates a square matrix with x on the main diagonal and 0’s on the off diagonal. The matrix of size N x N where N is length of the vector x. The command length(x) will return this length: >> length(x) Type: >> C=diag(x) The function diag can also be used to create matrices with entries above or below the main diagonal. If you want to create an N x N matrix with a vector z that is of size N-1 x 1, that is one above the main diagonal, you can type: >> z=[8; 4]; >> diag(z,1) to create a 3 x 3 matrix. Or if you want z to be 1 below the main diagonal type: >> diag(z,-1) For example, can you guess what you would get when you type in >> D=diag(x) + diag(x(1:2),1) - diag(x(2:3),-1) Now what we have a better understanding of how vectors/matrices work, we can proceed. Type >> clear 2

to delete all of the variables in our workspace. In Matlab, you work mostly with function files or script files (name.m). The difference is that in a function file, the first lines looks like: function [y,z]=name(x1,x2) where the outputs are y and z, the inputs are x1,x2. You can run this function in the command window by typing: >> y=name(x1,x2); or >> [y,z]=name(x1,x2); A script file is just a series of commands that you would type into the command line. There is a difference in ’scope’. The function file has a local ’scope’ (you can create new variables in there that will be deleted when you exit) versus in a script file, you have access to all the current variables in your workspace. Any changes to the variables you do in the script file will be changed in your workspace as well.

2. ODEs Matlab can solve large dimensional ODE systems with built-in solvers. You may want to consider using XPP if you are going to alter many parameters and do other dynamical system things (i.e., bifurcation diagrams, phaseplane analysis, etc.). Note that higher level programming is easier to implement in Matlab. Examples are complicated logical (if-else) statements in your system, weird boundary conditions, as well as look up tables [Sorry Bard]. You can use look-up tables in XPP as well, but I’m not sure how well it is documented. There are various methods (see Matlab’s documenation), but I will consider my favorite stiff solver: ode15s. For this workshop, we will assume the ODEs are first-order. For higher orders, we simply have to introduce auxiliary variables and enlarge the size of the system. Consider y ∈ R2 , and an ODE system of the form y˙1 = y1 · (a1 + a2 − y2 )

(1)

y˙2 = y1 − y2

(2)

with some initial conditions.

3

Open up the file rhs.m ; notice the syntax for this file. It basically specifies the right-hand side of the ODE system above. The parameters come after t (independent variable) and y (the functions y1 =y(1) and y2 =y(2). Fill in the appropriate commands for y˙1 =dy(1) and y˙2 =dy(2)..... ask for help if there is something unclear. Now open the file slv ode.m. Notice the initial conditions are random. Also the parameters for the time-length (t end) and a1 and a2 are already filled in (you can play with these later). Ignore the command options=..; you can explore this later by typing help ode15s in the command lind. You will write the syntax to solve the ODE [T,Y] = ode15s(@rhs,[0 t end],y0,options,a 1,a 2); Type the following in the command line: >> slv ode do this several times and see what happens. Before moving on to the next section, type: >> clear

3. PDEs We will discuss how to numerically solve partial differential equations of the form (x ∈ R1 and t ∈ [0, ∞)):

o o ∂u(x, t) ∂2 n ∂ n A(x, t)u(x, t) − 2 D(x)u(x, t) = 0. + ∂t ∂x ∂x

(3)

These are called ’advection-diffusion’ equations and arise all of the time in the physical and life sciences. This encompasses equations called ’Fokker-Planck’, ’Black-Scholes’, and a biological phenomena known as ’chemotaxis’ with some augmentations. For this workshop, we will consider the following advection-diffusion equation: o o ∂2 n ∂u(x, t) 1 ∂ n (µ(t) − x)u(x, t) − D 2 u(x, t) = 0 + ∂t τ ∂x ∂x

(4)

With zero boundary conditions: limx→±∞ u(x, t) = 0. o This can be understood intuitively n ∂ as follows. The advection term τ1 ∂x (µ(t) − x)u(x, t) causes the ’mass’ at x 6= µ(t) to

move (or advect) toward µ(t) with a velocity of τ1 . Thus, if τ is smaller, the ’mass’ moves

faster towards µ(t). Of course, the value µ(t) is changingnin timeo which adds to the ∂2 complications of these dynamics. The diffusion term D ∂x u(x, t) spreads or smears 2 4

out the mass equally in all directions (2 directions here since x ∈ R1 ) with a coefficient of

D that is constant here. Think of heat or gas diffusing on a line. Many different systems have a D that depends on x: D(x) (equation 3), which also complicates things.

3.1. Implementation Let us first rewrite the PDE o ∂u(x, t) ∂ 2 u(x, t) µ(t) ∂u(x, t) 1 ∂ n xu(x, t) + D = − ∂t τ ∂x ∂x2 τ ∂x

(5)

On a computer, we can’t go out to t = ∞, so we pick a large enough Tend = 50 say.

Discretize time with a mesh size ∆t help make matrix You should be able to see the form of the function and description. Now type >> [M0,Md,Am,x]=make matrix(25,200); Now we have our matrices in the workspace. Specify tau, D, and the end of time T end with the commands >> tau=2; >> D=5; >> T end=30; SOLVE THE PDE with the command: 7

[t,x,u,mu t]=solve advDif(tau,D,T end,M0,Md,Am,x); We will look at movies of our solutions in a bit.....

3.2. Monte Carlo Simulations Don’t be alarmed if you do not know any stochastic differential equations or statistical physics. But.. The advection-diffusion equation (4) describes how the probability density u(x, t) of a random variable X(t) evolves. The equation is: X˙ = µ(t) − X +

√

2Dξ(t)

where ξ(t) is a white noise process (i.e. a Gaussian random variables for all t with 0 mean and variance of 1; where each Gaussian random variable is uncorrelated in time). Since we took u(x, 0) = χ[0,1] , we can think of the initial condition X(0) as a random variable itself with uniform distribution. In the command line, type: >> avgX = monteCarlo sims(t,x,tau,D,mu t); After it is finished running, type: >> movie scpt to see a movie of how well the solution u(x, t) approximates the noisy system. The match is pretty good despite the crudeness of the methods (first order). Other issues to be aware of: i) mesh size isn’t large enough (in x or t, ii) conserving probability with weird boundary conditions, etc. You can go back and play with the function mu t in the file solve advDif if you are comfortable with the syntax, but you MUST enlarge your domain and mesh size ∆x appropriately!

4. Nonlinear PDEs Nonlinear PDEs are a mess to solve. Here is all I have to say about it. The instabilities that arise generally call for an implicit method (at the very least) to have any hope of a reasonable numerical solution. The following sections assume you want an implicit time-stepping solver.

8

4.1. Method of Lines If your PDE is nonlinear in u(x, t), then you have to use a fancier method than what has been outlined. One method is to write the PDE as a large system of ODEs and use the Method of Lines (Google it). This basically involves discretizing your domain x ∈ [−M, M] into N pieces of width ∆x = 2M/N, and letting u(xj , t) = u(−M + R −M +∆xj 1 ∆x(j − 1), t) represent ∆x u(x, t) dx. Note here that xj = −M + ∆x(j − 1) −M +∆x(j−1) for j = 1, 2, .., N. We still have

∂u ∂t

=

u(x,tk )−u(x,tk−1 ) . ∆t

The system of ODEs can be written down if you discretize the advection and diffusion parts with a finite difference method. For example, you could use a centered-difference method: o A(x , t)u(x , t) − A(x , t)u(x , t) ∂ n j+1 j+1 j−1 j−1 A(xj , t)u(xj , t) = ∂x 2∆x o D(x )u(x , t) − 2D(x , t)u(x ) + D(x )u(x , t) ∂2 n j+1 j+1 j j j−1 j−1 D(x )u(x , t) = j j ∂x2 ∆x2 Care must be taken with the boundary conditions. Again, you might want to use XPP if your boundary conditions are not nasty. There are many good nonlinear ODE solvers once you have properly discretized the system (e.g., CVODE, ode15s, etc.).

4.2. Newton’s Method Alternatively, you could use Newton’s Method to resolve the nonlinearity at each time step. This involves finding u(x, tk ) given that you know u(x, tk ). There are many deep and troubling issues such as uniqueness of roots, stability, etc. This is not unexpected with nonlinear PDEs.

9