CSCI 2227, Introduction to Scientific Computation
Prof. Alvarez



Partial Differential Equations (PDE)

	PDE are like ODE, in that unknown is a function

		ODE: seek u(t) such that u and the t-slope of u satisfy
		an equation of the form 

			du/dt(t) = f(t, u(t)), 
		
		and a pointwise condition 

			u(0) = u0

	PDE are unlike ODE, in that unknown is function of > 1 vars

		PDE: seek u(x,t) such that u, the t-slope of u, and the
		x-slope of u satisfy an equation, perhaps of the form 

			c du/dt(x,t) = f(x,t,u(x,t),du/dx(x,t))

		(higher-order derivatives may also enter)

		and boundary conditions of the form

			u(x,0) = f(x)  or  u(a,t) = d

	PDE provide continuous (distributed-parameter) models in n-D



Examples of PDE

	1) First-order convection equation

		du(x,t)/dt = c du(x,t)/dx

	subject to initial condition

		u(x,0) = f(x)

	2) Second-order diffusion equation

		du(x,t)/dt = s d^2u(x,t)/dx^2

		(s = diffusion constant)

	subject to initial condition

		u(x,0) = f(x)

	and boundary conditions

		u(a,t) = A,  u(b,t) = B

	3) Second-order electrostatic equation (diffusion equilibrium)

		d^2u(x,t)/dx^2 = 0

	subject to boundary conditions

		u(a,t) = A,  u(b,t) = B

	4) Black-Scholes option pricing

		dV(s,t)/dt + w^2 s^2/2 d^2V(s,t)ds^2 
		+ rs dV(s,t)/ds - rV(s,t) 
		= 0

		s = stock price
		t = time
		V(s,t) = fair price (value) of option
		w = volatility
		r = risk-free interest rate

	subject to boundary conditions such as (for call option):

		V(s,T) = max(s - K, 0),  
		(K is strike price at maturity (t=T))

		V(0,t) = 0 for all t



Finite-difference approximation of PDE solutions in 1 spatial dimension

	Based on finite-difference approximations of derivatives

		du(x,t)/dt approximated by 

			(u(x,t+dt) - u(x,t)) / dt

		du(x,t)/dx approximated by 

			(u(x+dx,t) - u(x,t)) / dx

		d^2u(x,t)/dx^2 approximated by 

			(u(x+dx,t) - 2u(x,t) + u(x-dx,t)) / dx^2

	For example, consider the diffusion equation:

		du(x,t)/dt = s^2 d^2u(x,t)/dx^2
		
	Finite difference approximation is:

		(u(x,t+dt) - u(x,t)) / dt 
		= s^2 (u(x+dx,t) - 2u(x,t) + u(x-dx,t)) / dx^2



Finite-differences as forward time-stepping schemes

	Can rewrite diffusion approximation in forward time as:

		u(x,t+dt) = u(x,t) 
			+ s^2 dt / dx^2 (u(x+dx,t) - 2u(x,t) + u(x-dx,t))

	(note that right-hand side only involves values at time t)

	Refer to grid points (i*dx, j*dt) by integer pairs (i,j)

	Define discrete version of u on integer pairs (i,j)

	Finite-difference version of diffusion PDE becomes:

		u(i,j+1) = u(i,j) 
			+ s^2 dt / dx^2 (u(i+1,j) - 2u(i,j) + u(i-1,j))



Finite-differences example in MATLAB

	We discussed the convection equation approximation in class

	using MATLAB function code [t x u] = convection(deltaX, deltaT)

		I posted a version on the CS227 homepage

	First run for deltaX = 0.1 = deltaT

	Results should be better for smaller deltaX, right?

	Run again with deltaX = 0.05, deltaT = 0.1

		What happens?



Stability of finite difference PDE approximation of convection PDE

	View u(x,t) as a function of x for each fixed t (time slice)

	Describe u(-,j) by its discrete Fourier transform for fixed j	

	Use the PDE to find DFT of u(-,j+1) in terms of DFT of u(-,j)

	Make use of the following important DFT shift identity:

		DFT( u(. + s) ) = e^isf2pi/N * DFT( u(.) )

	Convection equation yields:

		u^(f,j+1) = u^(f,j) [ 1 + dt/dx( e^i f 2pi/N - 1 ) ] 

	Modulus of factor on right (in square brackets) determines stability

	Note that factor traces out circle tangent to Im axis at 0

	Maximum modulus is either 

		1 (for f=0), if dt/dx < 1, 

		or

		2dt/dx - 1 (for f=N/2), if dt/dx >= 1

	Instability occurs only in the latter case, if dt/dx > 1

	Notice that maximum modulus occurs for frequency index N/2,
	which corresponds to a half-period of 1
 


Convection finite-difference example, revisited

	Examine the final time-slice in unstable case closely

	>> [t x u] = convection(.05,.0625);
	>> finalSlice = u(:,17);
	>> plot(finalSlice(180:190));

	Notice that oscillation has a half-period of 1, which is precisely 
	the mode of maximum DFT ratio modulus according to our analysis