Samir Hamdi
From Scholarpedia
Curator of ScholarpediaCurator Index: 0(Difference between revisions)
Samir Hamdi (Talk | contribs) m |
Samir Hamdi (Talk | contribs) m |
||
| Line 160: | Line 160: | ||
* Selected numerical results are stored for subsequent tabular and plotted output | * Selected numerical results are stored for subsequent tabular and plotted output | ||
<code><pre> | <code><pre> | ||
| − | |||
% Store numerical and analytical solutions, errors at x = 1/2 | % Store numerical and analytical solutions, errors at x = 1/2 | ||
n2=(n-1)/2.0+1; | n2=(n-1)/2.0+1; | ||
| Line 173: | Line 172: | ||
* Selected tabular numerical output is first displayed | * Selected tabular numerical output is first displayed | ||
<code><pre> | <code><pre> | ||
| − | |||
% Display selected output | % Display selected output | ||
fprintf('\n mf = %2d abstol = %8.1e reltol = %8.1e\n',... | fprintf('\n mf = %2d abstol = %8.1e reltol = %8.1e\n',... | ||
| Line 186: | Line 184: | ||
The output from this code is | The output from this code is | ||
| + | |||
| + | <code><pre> | ||
| + | mf = 1 abstol = 1.0e-004 reltol = 1.0e-004 | ||
| + | |||
| + | t u(0.5,t) u_anal(0.5,t) err u(0.5,t) | ||
| + | 0.000 0.707107 0.707107 0.0000000 | ||
| + | 0.625 0.151387 0.151268 0.0001182 | ||
| + | 1.250 0.032370 0.032360 0.0000093 | ||
| + | 1.875 0.006894 0.006923 -0.0000283 | ||
| + | 2.500 0.001472 0.001481 -0.0000091 | ||
| + | ncall = 85 | ||
| + | |||
| + | </pre></code> | ||
Revision as of 16:40, 25 September 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
%
% 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;
We can note the following points about this main program:
- After declaring some parameters global so that they can be shared with other routines called via this main program, initial condition (2) is computed over a 21-point grid in \( x \)
% 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
- The independent variable \( t \) is defined over the interval \( 0 \leq t\leq 2.5 \) again, a 21-point grid is used.
%
% Independent variable for ODE integration
t0=0.0;
tf=2.5;
tout=linspace(t0,tf,n);
nout=n;
ncall=0;
- The 21 ODEs are then integrated by a call to the Matlab integrator ode15s
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
Three cases are programmed corresponding to mf = 1, 2, 3 , for which three different ODEs routines, pde_1 , pde_2 , and pde_3 are called (these routines are discussed subsequently). The variable ndss refers to a library of differentiation routines for use in the MOL solution of PDEs; the use of ndss is illustrated in the subsequent discussion. Note that a stiff integrator, ode15s , was selected because the 21 ODEs are sufficiently stiff that a nonstiff integrator results in a large number of calls to the ODE routine.
- Selected numerical results are stored for subsequent tabular and plotted output
% 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
- Selected tabular numerical output is first displayed
% 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);
The output from this code is
mf = 1 abstol = 1.0e-004 reltol = 1.0e-004
t u(0.5,t) u_anal(0.5,t) err u(0.5,t)
0.000 0.707107 0.707107 0.0000000
0.625 0.151387 0.151268 0.0001182
1.250 0.032370 0.032360 0.0000093
1.875 0.006894 0.006923 -0.0000283
2.500 0.001472 0.001481 -0.0000091
ncall = 85
