pde_notes04

# pde_notes04 - u x 0 = φ x = ± 1 | x |/a if | x |< a...

This preview shows page 1. Sign up to view the full content.

Notes for TUT course MAT-51316 Robert Pich´e 10.9.2007 4 Plot Solution of 1D Diffusion IVP The solution of the one dimensional diffusion equation u t = ku xx with initial condition u ( x, 0) = φ ( x ) is given by the convolution integral u ( x, t ) = Z -∞ S ( x - y, t ) φ ( y ) dy where S ( x ) = 1 4 πkt e - x 2 / (4 kt ) . To compute the convolution integral, ﬁrst choose equally spaced points y 1: N such that the integrand is negligible outside the interval [ y 1 , y N ] . Then, applying the midpoint formula Z b a F ( y ) dy ( b - a ) F ( a + b 2 ) . to approximate the convolution integral in each subinterval gives u ( x, t ) Z y N y 1 S ( x - y, t ) φ ( y ) dy N X j =1 hS ( x - y j - h 2 , t ) φ ( y j + h 2 ) where h = ( y N - y 1 ) / ( N - 1) . The following code computes and plots the solution for the initial condition
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: u ( x, 0) = φ ( x ) = ± 1- | x | /a if | x | < a otherwise a=1; b=1; k=1; % scaling for x, u, t [email protected](x,a,b) b*(1-abs(x)/a).*(abs(x)<a); [email protected](x,t,k) exp(-x.ˆ2 ./(4*k*t))./sqrt(4*pi*k*t); nx=41; nt=20; ny=11; x=linspace(-4*a,4*a,nx); y=linspace(-a,a,ny); t=linspace(0,3*aˆ2/k,nt); [xx,yy]=meshgrid(x,y); h=y(2)-y(1); f=phi(y+h/2,a,b); u=zeros(nt,nx); u(1,:)=phi(x,a,b); for it=2:nt u(it,:)=. .. h*f*S(xx-yy-h/2,t(it),k); end surf(x,t,u) − 4 − 2 2 4 0.5 1 1.5 2 2.5 3 0.2 0.4 0.6 0.8 1 1...
View Full Document

{[ snackBarMessage ]}

Ask a homework question - tutors are online