|
|
| Line 1: |
Line 1: |
| − | 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
| |
| − |
| |
| − | <code><pre>
| |
| − |
| |
| − | %
| |
| − | % 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;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − | We can note the following points about this main program:
| |
| − |
| |
| − | * After declaring some parameters <tt> global </tt> 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 <math> x </math>
| |
| − |
| |
| − | <code><pre>
| |
| − | % 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
| |
| − | </pre></code>
| |
| − |
| |
| − | * The independent variable <math> t </math> is defined over the interval <math> 0 \leq t\leq 2.5 </math> again, a 21-point grid is used.
| |
| − | <code><pre>
| |
| − | %
| |
| − | % Independent variable for ODE integration
| |
| − | t0=0.0;
| |
| − | tf=2.5;
| |
| − | tout=linspace(t0,tf,n);
| |
| − | nout=n;
| |
| − | ncall=0;
| |
| − | </pre></code>
| |
| − |
| |
| − | * The 21 ODEs are then integrated by a call to the Matlab integrator <tt> ode15s </tt>
| |
| − | <code><pre>
| |
| − | 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
| |
| − | </pre></code>
| |
| − | Three cases are programmed corresponding to <tt> mf = 1, 2, 3 </tt> , for
| |
| − | which three different ODEs routines, <tt> pde_1 </tt>, <tt> pde_2 </tt>, and
| |
| − | <tt> pde_3 </tt> are called (these routines are discussed subsequently). The
| |
| − | variable <tt> ndss </tt> refers to a library of differentiation routines
| |
| − | for use in the MOL solution of PDEs; the use of <tt> ndss </tt> is illustrated
| |
| − | in the subsequent discussion. Note that a stiff integrator, <tt> ode15s </tt>,
| |
| − | 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
| |
| − | <code><pre>
| |
| − | % 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
| |
| − | </pre></code>
| |
| − |
| |
| − | * Selected tabular numerical output is first displayed
| |
| − | <code><pre>
| |
| − | % 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);
| |
| − | </pre></code>
| |
| − |
| |
| − | 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>
| |
| − |
| |
| − | This output indicates that the MOL solution agrees with the analytical
| |
| − | solution to at least three significant figures. Also, <tt> ode15s </tt>
| |
| − | calls the derivative routine only 85 times (in contrast with the nonstiff
| |
| − | integrator <tt> ode45 </tt> 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
| |
| − |
| |
| − | <code><pre>
| |
| − | % 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
| |
| − | </pre></code>
| |
| − |
| |
| − |
| |
| − | The plotted error output below indicates that the error in the MOL
| |
| − | solution varied between approximately <math> -3\times 10^{-5} </math> and
| |
| − | <math>16\times 10^{-5}</math> which is not quite within the error range specified in the
| |
| − | program
| |
| − |
| |
| − | <code><pre>
| |
| − |
| |
| − | reltol=1.0e-04; abstol=1.0e-04;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − | The fact that these error tolerances were not satisfied does not necessarily
| |
| − | mean that <tt> ode15s </tt> failed to adjust the integration interval to meet
| |
| − | these error tolerances. Rather, the error of approximately <math> 1.6 \times 10^{-4} </math>
| |
| − | is due to the limited accuracy of the second order FD approximation of <math> \frac{\partial ^2 u}
| |
| − | {\partial x^2} </math> programmed in <tt> pde_1 </tt>. This conclusion is confirmed when
| |
| − | the main program calls <tt> pde_2 </tt> (for <tt> mf = 2 </tt> ) or <tt> pde_3 </tt>
| |
| − | (for <tt> mf = 3 </tt>) as discussed subsequently; these two routines have FD approximations that are more
| |
| − | accurate than in <tt> pde_1 </tt> 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 <math> t </math>
| |
| − | (by <tt> ode15s </tt>) and (b)
| |
| − | errors due to the approximation of the spatial derivatives such as
| |
| − | <math> \frac{\partial ^2 u}{\partial x^2} </math> programmed in the derivative routine such
| |
| − | as <tt> pde_1 </tt>. 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 <math> x </math> were not sufficient when using the second order FDs
| |
| − | in <tt> pde_1 </tt>.
| |
| − |
| |
| − | 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
| |
| − |
| |
| − | <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;
| |
| − | </pre></code>
| |
| − |
| |
| − | We can note the following points about <tt> tt pde_1 </tt>:
| |
| − |
| |
| − | * Some problem parameters are first defined
| |
| − |
| |
| − | <code><pre>
| |
| − |
| |
| − | function ut=pde_1(t,u)
| |
| − | %
| |
| − | % Problem parameters
| |
| − | global ncall
| |
| − | xl=0.0;
| |
| − | xu=1.0;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − | <tt> xl </tt> and <tt> xu </tt> could have also been set in the main program and passed
| |
| − | to <tt> pde_1 </tt> as global variables. The defining statement at the beginning
| |
| − | of <tt> pde_1 </tt> indicates that the independent variable <tt> t </tt> and dependent
| |
| − | variable vector <tt> u </tt> are inputs to <tt> pde_1 </tt>, while the output is the vector
| |
| − | of <tt> tt t </tt> dervatives, <tt> ut </tt>; in other words, all of the <tt> n </tt> ODE derivatives
| |
| − | in <tt> t </tt> must be defined in <tt> pde_1 </tt>.
| |
| − |
| |
| − | * The finite difference approximation of eq. (1) is then programmed
| |
| − |
| |
| − | <code><pre>
| |
| − | % 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';
| |
| − | </pre></code>
| |
| − |
| |
| − | The number of ODEs (21) is determined by the length command <tt> n=length(u); </tt>
| |
| − | so that the programming is general (the number of ODEs can easily be
| |
| − | changed in the main program). The square of the FD interval, <tt> dx2 </tt>,
| |
| − | is then computed.
| |
| − |
| |
| − | * The MOL programming of the 21 ODEs is done in the <tt> for </tt> loop. For BC
| |
| − | (3), the coding is
| |
| − |
| |
| − | <code><pre>
| |
| − |
| |
| − | if(i==1) ut(i)=0.0;
| |
| − | </pre></code>
| |
| − |
| |
| − | since the value of <math> u(x=0,t) = 0 </math> does not change after being set as
| |
| − | an initial condition in the main program (and therefore its time derivative
| |
| − | is zero).
| |
| − |
| |
| − | * For BC (4), the coding is
| |
| − |
| |
| − | <code><pre>
| |
| − |
| |
| − | elseif(i==n) ut(i)=2.0*(u(i-1)-u(i))/dx2;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − | which follows directly from the FD approximation of BC (4) (eq. (37))
| |
| − |
| |
| − | :<math>
| |
| − | u_x \approx \frac{u(i+1) - u(i-1)}{\Delta x} = 0
| |
| − | </math>
| |
| − |
| |
| − | or with <math> i = n </math>
| |
| − |
| |
| − | :<math>
| |
| − | u(n+1) = u(n-1)
| |
| − | </math>
| |
| − |
| |
| − | Note that the ''fictitious value'' <math> u(n+1) </math> can then be replaced in the ODE
| |
| − | at <math>i=n</math> by <math>u(n-1) </math>.
| |
| − |
| |
| − | * For the remaining interior points, the programming is
| |
| − |
| |
| − | <code><pre>
| |
| − | else ut(i)=(u(i+1)-2.0*u(i)+u(i-1))/dx2;
| |
| − | </pre></code>
| |
| − |
| |
| − |
| |
| − |
| |
| − | which follows from the FD approximation of the second derivative (eq. (33))
| |
| − |
| |
| − |
| |
| − | :<math>
| |
| − | u_{xx} \approx \frac{(u(i+1) - 2u(i) + u(i-1))}{\Delta x^2}
| |
| − | </math>
| |
| − |
| |
| − |
| |
| − | * Since the Matlab ODE integrators require a column vector of derivatives, a final transpose of <tt> ut </tt> is required
| |
| − |
| |
| − | <code><pre>
| |
| − |
| |
| − | ut=ut';
| |
| − | %
| |
| − | % Increment calls to pde_1
| |
| − | ncall=ncall+1;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − |
| |
| − | Finally, the number of calls to <tt> pde_1 </tt> is incremented so that at the
| |
| − | end of the solution, the value of <tt> ncall </tt> displayed by the main program
| |
| − | gives an indication of the computational effort required to produce
| |
| − | the entire solution. The numerical and graphical output for this case
| |
| − | (<tt>mf=1</tt>) was discussed subsequently.
| |
| − |
| |
| − | For <tt>mf=2</tt>, function <tt> pde_2 </tt> is called by <tt> ode15s </tt>
| |
| − |
| |
| − | <code><pre>
| |
| − |
| |
| − | function ut=pde_2(t,u)
| |
| − | %
| |
| − | % Problem parameters
| |
| − | global ncall ndss
| |
| − | xl=0.0;
| |
| − | xu=1.0;
| |
| − | %
| |
| − | % BC at x = 0 (Dirichlet)
| |
| − | u(1)=0.0;
| |
| − | %
| |
| − | % Calculate ux
| |
| − | n=length(u);
| |
| − | if (ndss== 2) ux=dss002(xl,xu,n,u); % second order
| |
| − | elseif(ndss== 4) ux=dss004(xl,xu,n,u); % fourth order
| |
| − | elseif(ndss== 6) ux=dss006(xl,xu,n,u); % sixth order
| |
| − | elseif(ndss== 8) ux=dss008(xl,xu,n,u); % eighth order
| |
| − | elseif(ndss==10) ux=dss010(xl,xu,n,u); % tenth order
| |
| − | end
| |
| − | %
| |
| − | % BC at x = 1 (Neumann)
| |
| − | ux(n)=0.0;
| |
| − | %
| |
| − | % Calculate uxx
| |
| − | if (ndss== 2) uxx=dss002(xl,xu,n,ux); % second order
| |
| − | elseif(ndss== 4) uxx=dss004(xl,xu,n,ux); % fourth order
| |
| − | elseif(ndss== 6) uxx=dss006(xl,xu,n,ux); % sixth order
| |
| − | elseif(ndss== 8) uxx=dss008(xl,xu,n,ux); % eighth order
| |
| − | elseif(ndss==10) uxx=dss010(xl,xu,n,ux); % tenth order
| |
| − | end
| |
| − | %
| |
| − | % PDE
| |
| − | ut=uxx';
| |
| − | ut(1)=0.0;
| |
| − | %
| |
| − | % Increment calls to pde_2
| |
| − | ncall=ncall+1;
| |
| − |
| |
| − | </pre></code>
| |
| − | We can note the following points about <tt> pde_2</tt>:
| |
| − |
| |
| − | * The initial statements are the same as in <tt> pde_1 </tt>. Then the Dirichlet
| |
| − | BC at <math> x = 0 </math> is programmed
| |
| − |
| |
| − | <code><pre>
| |
| − | % BC at x = 0 (Dirichlet)
| |
| − | u(1)=0.0;
| |
| − | </pre></code>
| |
| − |
| |
| − | Acually, the statement <tt> u(1)=0.0;</tt> has no effect since the dependent
| |
| − | variables can only be changed through their derivatives, i.e., <tt> ut(1) </tt>,
| |
| − | in the ODE derivative routine. This code was included just to serve
| |
| − | as a reminder of the BC at <math>x = 0<math>, which is programmed subsequently
| |
| − |
| |
| − | * The first order spatial derivative <math> \frac{\partial u}{\partial x}=u_{x} </math>
| |
| − | is then computed
| |
| − |
| |
| − | <code><pre>
| |
| − | % Calculate ux
| |
| − | n=length(u);
| |
| − | if (ndss==2) ux=dss002(xl,xu,n,u); % second order
| |
| − | elseif(ndss== 4) ux=dss004(xl,xu,n,u); % fourth order
| |
| − | elseif(ndss== 6) ux=dss006(xl,xu,n,u); % sixth order
| |
| − | elseif(ndss== 8) ux=dss008(xl,xu,n,u); % eighth order
| |
| − | elseif(ndss==10) ux=dss010(xl,xu,n,u); % tenth order
| |
| − | end
| |
| − | </pre></code>
| |
| − |
| |
| − |
| |
| − | Five library routines, <tt> dss002</tt> to <tt> dss010</tt>, are programmed that use second
| |
| − | order to tenth order FD approximations, respectively. Since <tt> ndss=4</tt>
| |
| − | is specified in the main program, <tt> dss004</tt> is used in the calculation
| |
| − | of <tt> ux</tt>.
| |
| − |
| |
| − | * BC (4) is then applied followed by the calculation of the second
| |
| − | order spatial derivative from the first order spatial derivative
| |
| − |
| |
| − | <code><pre>
| |
| − | % BC at x = 1 (Neumann)
| |
| − | ux(n)=0.0;
| |
| − | %
| |
| − | % Calculate uxx
| |
| − | if (ndss== 2) uxx=dss002(xl,xu,n,ux); % second order
| |
| − | elseif(ndss== 4) uxx=dss004(xl,xu,n,ux); % fourth order
| |
| − | elseif(ndss== 6) uxx=dss006(xl,xu,n,ux); % sixth order
| |
| − | elseif(ndss== 8) uxx=dss008(xl,xu,n,ux); % eighth order
| |
| − | elseif(ndss==10) uxx=dss010(xl,xu,n,ux); % tenth order
| |
| − | end
| |
| − | </pre></code>
| |
| − |
| |
| − | Again, <tt> dss004</tt> is called which is the usual procedure (the order of
| |
| − | the FD approximation generally is not changed in computing higher
| |
| − | order derivatives from lower order derivative, a process termed ''stagewise
| |
| − | differentiation'').
| |
| − |
| |
| − | * Finally, eq. (1) is programmed and the Dirichlet BC at <math> x = 0 </math>
| |
| − | (eq.(3)) is applied
| |
| − |
| |
| − | <code><pre>
| |
| − | % PDE
| |
| − | ut=uxx';
| |
| − | ut(1)=0.0;
| |
| − | %
| |
| − | % Increment calls to pde_2
| |
| − | ncall=ncall+1;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − | Note the similarity of the code to the PDE (eq. (1)), and also the
| |
| − | transpose required by <tt> ode15s</tt>.
| |
| − |
| |
| − | The numerical output for this case (<tt>mf=2</tt>) is
| |
| − |
| |
| − | <code><pre>
| |
| − | mf = 2 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.151267 0.151268 -0.0000013
| |
| − | 1.250 0.032318 0.032360 -0.0000418
| |
| − | 1.875 0.006878 0.006923 -0.0000446
| |
| − | 2.500 0.001467 0.001481 -0.0000138
| |
| − |
| |
| − | ncall = 62
| |
| − | </pre></code>
| |
| − |
| |
| − |
| |
| − | The plotted error output below indicates that the error in the MOL
| |
| − | solution varied between approximately <math>-5\times 10^{-5}</math> and
| |
| − | <math>3.2\times 10^{-5}</math> which is within the error range specified in the
| |
| − | program
| |
| − |
| |
| − | <code><pre>
| |
| − |
| |
| − | reltol=1.0e-04; abstol=1.0e-04;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − |
| |
| − | Thus, switching from the second order FDs in <tt>pde_1</tt> to fourth order
| |
| − | finite differences in <tt> pde_2</tt> reduced the ''spatial truncation error''
| |
| − | so that the MOL solution met the specified error tolerances.
| |
| − |
| |
| − | For <tt> mf = 3</tt>, function <tt> pde_3</tt> is called by <tt> ode15s </tt>
| |
| − |
| |
| − | <code><pre>
| |
| − | function ut=pde_3(t,u)
| |
| − | %
| |
| − | % Problem parameters
| |
| − | global ncall ndss
| |
| − | xl=0.0;
| |
| − | xu=1.0;
| |
| − | %
| |
| − | % BC at x = 0
| |
| − | u(1)=0.0;
| |
| − | %
| |
| − | % BC at x = 1
| |
| − | n=length(u);
| |
| − | ux(n)=0.0;
| |
| − | %
| |
| − | % Calculate uxx
| |
| − | nl=1; % Dirichlet
| |
| − | nu=2; % Neumann
| |
| − | if (ndss==42) uxx=dss042(xl,xu,n,u,ux,nl,nu); % second order
| |
| − | elseif(ndss==44) uxx=dss044(xl,xu,n,u,ux,nl,nu); % fourth order
| |
| − | elseif(ndss==46) uxx=dss046(xl,xu,n,u,ux,nl,nu); % sixth order
| |
| − | elseif(ndss==48) uxx=dss048(xl,xu,n,u,ux,nl,nu); % eighth order
| |
| − | elseif(ndss==50) uxx=dss050(xl,xu,n,u,ux,nl,nu); % tenth order
| |
| − | end
| |
| − | %
| |
| − | % PDE
| |
| − | ut=uxx';
| |
| − | ut(1)=0.0;
| |
| − | %
| |
| − | % Increment calls to pde_3
| |
| − | ncall=ncall+1;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − | We can note the following points about <tt>pde_3</tt>:
| |
| − |
| |
| − | * The initial statements are the same as in <tt> pde_1</tt>. Then the Dirichlet
| |
| − | BC at <math>x = 0</math> and the Neumann BC at <math>x = 1 </math> are programmed
| |
| − |
| |
| − | <code><pre>
| |
| − | function ut=pde_3(t,u)
| |
| − | %
| |
| − | % Problem parameters
| |
| − | global ncall ndss
| |
| − | xl=0.0;
| |
| − | xu=1.0;
| |
| − | %
| |
| − | % BC at x = 0
| |
| − | u(1)=0.0;
| |
| − | %
| |
| − | % BC at x = 1
| |
| − | n=length(u);
| |
| − | ux(n)=0.0;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − | Again, the statement <tt> u(1)=0.0 </tt>; has no effect (since the dependent
| |
| − | variables can only be changed through their derivatives, i.e., <tt> ut(1)</tt>,
| |
| − | in the ODE derivative routine). This code was included just to serve
| |
| − | as a reminder of the BC at <math> x = 0 </math>, which is programmed subsequently.
| |
| − |
| |
| − | * The second order spatial derivative <math>\frac{\partial^{2}u}{\partial x^{2}}=u_{xx}</math>
| |
| − | is then computed
| |
| − |
| |
| − | <code><pre>
| |
| − | % Calculate uxx
| |
| − | nl=1; % Dirichlet
| |
| − | nu=2; % Neumann
| |
| − | if (ndss==42) uxx=dss042(xl,xu,n,u,ux,nl,nu); % second order
| |
| − | elseif(ndss==44) uxx=dss044(xl,xu,n,u,ux,nl,nu); % fourth order
| |
| − | elseif(ndss==46) uxx=dss046(xl,xu,n,u,ux,nl,nu); % sixth order
| |
| − | elseif(ndss==48) uxx=dss048(xl,xu,n,u,ux,nl,nu); % eighth order
| |
| − | elseif(ndss==50) uxx=dss050(xl,xu,n,u,ux,nl,nu); % tenth order
| |
| − | end
| |
| − | </pre></code>
| |
| − |
| |
| − | Five library routines, <tt>dss042</tt> to <tt>dss050</tt>, are programmed that use second
| |
| − | order to tenth order FD approximations, respectively for a second
| |
| − | derivative. Since <tt> ndss=44</tt> is specified in the main program, <tt>dss0044</tt>
| |
| − | is used in the calculation of <tt> uxx</tt>. Also, these differentiation routines
| |
| − | have two parameters that specify the type of BCs: (a) <tt> nl = 1 </tt> or <tt> 2 </tt>
| |
| − | specify a Dirichlet or a Neumann BC, respectively, at the lower boundary
| |
| − | value of <math> x = xl (= 0)</math>; in this case, BC (3) is Dirichlet, so <tt> nl = 1</tt>
| |
| − | and (b) <tt> nu = 1</tt> or <tt> 2</tt> specify a Dirichlet or a Neumann BC, respectively,
| |
| − | at the upper boundary value of <math>x = xu (= 1)</math>; in this case, BC (4)
| |
| − | is Neumann, so <tt> nu = 2</tt>.
| |
| − |
| |
| − | * Finally, eq. (1) is programmed and the Dirichlet BC at <math>x = 0</math> (eq. (3)) is applied
| |
| − |
| |
| − | <code><pre>
| |
| − | %
| |
| − | % PDE
| |
| − | ut=uxx';
| |
| − | ut(1)=0.0;
| |
| − | %
| |
| − | % Increment calls to pde_3
| |
| − | ncall=ncall+1;
| |
| − |
| |
| − | </pre></code>
| |
| − |
| |
| − | Again, the transpose is required by <tt> ode15s</tt>.
| |
| − |
| |
| − | The numerical output for this case (<tt>mf=3</tt>) is
| |
| − |
| |
| − | <code><pre>
| |
| − |
| |
| − | mf = 3 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.151267 0.151268 -0.0000017
| |
| − | 1.250 0.032318 0.032360 -0.0000420
| |
| − | 1.875 0.006878 0.006923 -0.0000447
| |
| − | 2.500 0.001467 0.001481 -0.0000138
| |
| − |
| |
| − | ncall = 62
| |
| − | </pre></code>
| |
| − |
| |
| − | The plotted error output below indicates that the error in the MOL
| |
| − | solution varied between approximately <math>-4.8\times 10^{-5}</math> and
| |
| − | <math>3.2\times 10^{-5}</math> which is within the error range specified in the
| |
| − | program
| |
| − |
| |
| − | <code><pre>
| |
| − | reltol=1.0e-04; abstol=1.0e-04;
| |
| − | </pre></code>
| |
| − |
| |
| − | We conclude this example with the following observation:
| |
| − | As the solution approaches steady state, <math>t\rightarrow\infty</math>,<math>u_{t}\rightarrow 0</math> and
| |
| − | from eq. (1), <math>u_{xx}\rightarrow 0</math>. As the second derivative vanishes, the solution becomes
| |
| − |
| |
| − | <code><pre>
| |
| − | u_{xx}=0
| |
| − | \end{equation*}
| |
| − | \begin{equation*}
| |
| − | u_{x}=c_{1}
| |
| − | \end{equation*}
| |
| − | \begin{equation*}
| |
| − | u=c_{1}x+c_{2}
| |
| − | </pre></code>
| |
| − |
| |
| − | Thus, the steady state solution is linear in <math>x</math>, which can serve as
| |
| − | another check on the numerical solution (for BCs (3) and (4), <math>c_{1}=c_{2}=0</math>
| |
| − | and thus at steady state <math>u=0,</math> which also follows from the analytical
| |
| − | solution, eq. (42)). This type of special case analysis is often useful
| |
| − | in checking a numerical solution. In addition to mathematical conditions
| |
| − | such as the linear dependency on <math>x</math>, physical conditions can frequently
| |
| − | also be used to check solutions, e.g., conservation of mass, momentum
| |
| − | and energy.
| |
| − |
| |
| − | Through this example application we have attempted to illustrate
| |
| − | the basic steps of MOL/PDE analysis to arrive at a numerical solution
| |
| − | of acceptable accuracy. We have also presented some basic ideas for
| |
| − | assessing accuracy with respect to time and space (e.g., <math>t</math> and <math>x</math>).
| |
| − | More advanced applications (e.g., problems expressed as systems of
| |
| − | nonlinear PDEs) can be analyzed using the same ideas discussed in this
| |
| − | example. The authors will be glad to assist by providing codes for
| |
| − | other 1D, 2D and 3D PDE problems.
| |
