pde_notes04

pde_notes04 - u ( x, 0) = ( x ) = 1- | x | /a if | x |...

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

View Full Document Right Arrow Icon
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, first 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
Background image of page 1
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 phi=@(x,a,b) b*(1-abs(x)/a).*(abs(x)<a); S=@(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*a2/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

Ask a homework question - tutors are online