Samir Hamdi
From Scholarpedia
Curator of ScholarpediaCurator Index: 0(Difference between revisions)
Samir Hamdi (Talk | contribs) m |
Samir Hamdi (Talk | contribs) m |
||
| Line 137: | Line 137: | ||
</pre></code> | </pre></code> | ||
| − | * The | + | * The 21 ODEs are then integrated by a call to the Matlab integrator <tt> ode15s </tt></math> |
Revision as of 14:42, 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 </math>
