Samir Hamdi

Revision as of 18:48, 23 August 2007

We now revisit the PDE problem of eqs. (1) to (4) to illustrate the details of constructing a MOL code and to discuss the numerical and graphical output from the code.

The specific initial condition of eq. (2) will be taken as:

<math eq_2>

u(x,t=0) = sin(\pi x/2) </math>

The analytical solution to eqs. (1) to (4) which we use to evaluate the numerical MOL output is

<math eq_2>

u(x,t)=e^{-(\pi^{2}/4)t}sin(\pi x/2) </math>

A main program in Matlab for the MOL solution of eqs. (1) to (4) with the analytical solution, eq. (42), included for comparison with the MOL solution, is listed below

\begin{verbatim} % % Clear previous files clear all clc % % Parameters shared with the ODE routine global ncall ndss % % Initial condition n=21; for i=1:n u0(i)=sin((pi/2.0)*(i-1)/(n-1)); end % % Independent variable for ODE integration t0=0.0; tf=2.5; tout=linspace(t0,tf,n); nout=n; ncall=0; % % ODE integration mf=1; reltol=1.0e-04; abstol=1.0e-04; options=odeset('RelTol',reltol,'AbsTol',abstol); if(mf==1) % explicit FDs [t,u]=ode15s(@pde_1,tout,u0,options); end if(mf==2) ndss=4; % ndss = 2, 4, 6, 8 or 10 required [t,u]=ode15s(@pde_2,tout,u0,options); end if(mf==3) ndss=44; % ndss = 42, 44, 46, 48 or 50 required [t,u]=ode15s(@pde_3,tout,u0,options); end % % Store numerical and analytical solutions, errors at x = 1/2 n2=(n-1)/2.0+1; sine=sin(pi/2.0*0.5); for i=1:nout u_plot(i)=u(i,n2); u_anal(i)=exp(-pi^2/4.0*t(i))*sine; err_plot(i)=u_plot(i)-u_anal(i); end % % Display selected output fprintf('\n mf = %2d abstol = %8.1e reltol = %8.1e\n',... mf,abstol,reltol); fprintf('\n t u(0.5,t) u_anal(0.5,t) err u(0.5,t)\n'); for i=1:5:nout fprintf('%6.3f%15.6f%15.6f%15.7f\n',... t(i),u_plot(i),u_anal(i),err_plot(i)); end fprintf('\n ncall = %4d\n',ncall); % % Plot numerical solution and errors at x = 1/2 figure(1); subplot(1,2,1) plot(t,u_plot); axis tight title('u(0.5,t) vs t'); xlabel('t'); ylabel('u(0.5,t)') subplot(1,2,2) plot(t,err_plot); axis tight title('Err u(0.5,t) vs t'); xlabel('t'); ylabel('Err u(0.5,t)') print -deps pde.eps; print -dps pde.ps % % Plot numerical solution in 3D perspective figure(2); colormap('Gray'); C=ones(n); g=linspace(0,1,n); % For distance x h1 = waterfall(t,g,u',C); axis('tight'); grid off xlabel('t, time') ylabel('x, distance') zlabel('u(x,t)') s1 = sprintf('Diffusion Equation - MOL Solution'); sTmp = sprintf('u(x,0) = sin(\\pi x/2 )'); s2 = sprintf('Initial condition: %s', sTmp); title([{s1}, {s2}], 'fontsize', 12); v = [0.8616 -0.5076 0.0000 -0.1770 0.3712 0.6301 0.6820 -0.8417 0.3462 0.5876 -0.7313 8.5590 0 0 0 1.0000]; view(v); rotate3d on; \end{verbatim}