Samir Hamdi
Samir Hamdi (Talk | contribs) m |
Samir Hamdi (Talk | contribs) m |
||
| Line 261: | Line 261: | ||
* A 3D plot is also produced | * A 3D plot is also produced | ||
| + | |||
| + | <code><pre> | ||
| + | % 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; | ||
| + | </pre></code> | ||
| + | |||
| + | The plotted output below clearly indicates the origin of the ''lines'' in the ''method of lines'' | ||
| + | |||
| + | The programming of the approximating MOL/ODEs is in one of the three | ||
| + | routines called by <tt> ode15s </tt>. We now consider each of these routines. | ||
| + | For <tt> mf = 1, pde_1 </tt> calls function <tt> ut=pde_1(t,u) <tt> | ||
| + | |||
| + | <code><pre> | ||
| + | |||
| + | function ut=pde_1(t,u) | ||
| + | % | ||
| + | % Problem parameters | ||
| + | global ncall | ||
| + | xl=0.0; | ||
| + | xu=1.0; | ||
| + | % | ||
| + | % PDE | ||
| + | n=length(u); | ||
| + | dx2=((xu-xl)/(n-1))^2; | ||
| + | for i=1:n | ||
| + | if(i==1) ut(i)=0.0; | ||
| + | elseif(i==n) ut(i)=2.0*(u(i-1)-u(i))/dx2; | ||
| + | else ut(i)=(u(i+1)-2.0*u(i)+u(i-1))/dx2; | ||
| + | end | ||
| + | end | ||
| + | ut=ut'; | ||
| + | % | ||
| + | % Increment calls to pde_1 | ||
| + | ncall=ncall+1; | ||
| + | <code><pre></pre> | ||
Revision as of 17:33, 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
This output indicates that the MOL solution agrees with the analytical solution to at least three significant figures. Also, ode15s calls the derivative routine only 85 times (in contrast with the nonstiff integrator ode45 that requires approximately 5000 - 10000 calls, which clearly indicates the advantage of a stiff integrator for this problem).
- The MOL solution and its error (computed from the analytical solution) are plotted
% 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
The plotted error output below indicates that the error in the MOL
solution varied between approximately \( -3\times 10^{-5} \) and
\(16\times 10^{-5}\) which is not quite within the error range specified in the
program
reltol=1.0e-04; abstol=1.0e-04;
The fact that these error tolerances were not satisfied does not necessarily mean that ode15s failed to adjust the integration interval to meet these error tolerances. Rather, the error of approximately \( 1.6 \times 10^{-4} \) is due to the limited accuracy of the second order FD approximation of \( \frac{\partial ^2 u} {\partial x^2} \) programmed in pde_1 . This conclusion is confirmed when the main program calls pde_2 (for mf = 2 ) or pde_3 (for mf = 3 ) as discussed subsequently; these two routines have FD approximations that are more accurate than in pde_1 so the errors fall below the specified tolerances.
This analysis indicates that two sources of errors result from the MOL solution
of PDEs such as eq. (1): (a) errors due to the integration in \( t \)
(by ode15s ) and (b)
errors due to the approximation of the spatial derivatives such as
\( \frac{\partial ^2 u}{\partial x^2} \) programmed in the derivative routine such
as pde_1 . In other words, we have to be attentive to integration
errors in the initial and boundary value independent variables.
In summary, a comparison of the numerical and analytical solutions indicates that 21
grid points in \( x \) were not sufficient when using the second order FDs
in pde_1 .
However, in general, we will not have an analytical solution such as eq. (42) to determine if the number of spatial grid points is adequate. In this case, some experimentation with the number of grid points, and the observation of the resulting solutions to infer the degree of accuracy or spatial convergence, may be required.
- A 3D plot is also produced
% 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;
The plotted output below clearly indicates the origin of the lines in the method of lines
The programming of the approximating MOL/ODEs is in one of the three routines called by ode15s . We now consider each of these routines. For mf = 1, pde_1 calls function ut=pde_1(t,u)
function ut=pde_1(t,u)
%
% Problem parameters
global ncall
xl=0.0;
xu=1.0;
%
% PDE
n=length(u);
dx2=((xu-xl)/(n-1))^2;
for i=1:n
if(i==1) ut(i)=0.0;
elseif(i==n) ut(i)=2.0*(u(i-1)-u(i))/dx2;
else ut(i)=(u(i+1)-2.0*u(i)+u(i-1))/dx2;
end
end
ut=ut';
%
% Increment calls to pde_1
ncall=ncall+1;
<code><pre>
