Partial differential equation/Approximate and Numerical Methods/Appendices
Contents |
Approximate and Numerical Methods: Appendices
Appendix 1
The MATLAB program that produced the results described in Section 1 for parabolic PDE () follows.
clear all; clc
D=1; mu=1; A=1; B=1;
nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11;
nt=50; nout=10; h=0.001; t=0;
for i=1:nx; w(i)=A*exp(1/D*(mu^2*t-mu*x(i)))+B; end
for i1=1:nout
for i2=1:nt
for i=1:nx
if(i==1)w(1)=A*exp(1/D*(mu^2*t-mu*x(1)))+B; wt(1)=0;
elseif(i==nx)wf=w(nx-1)+(2*dx)*A*(-mu/D)*exp(1/D*(mu^2*t-mu*x(nx)));
wt(nx)=D*(wf-2*w(nx)+w(nx-1))/dx^2;
else wt(i)=D*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
end
w=w+wt*h; t=t+h;
end
fprintf('t = %7.4f w(x=xl/2,t) = %7.4f wa(x=xl/2,t) = %7.4f\n', ...
t, w(nx2), A*exp(1/D*(mu^2*t-mu*x(nx2)))+B);
end
Listing 1: MATLAB program for the numerical solution of Eqs (), (), (), and ()
We can note the following points about this program:
- Any previous files are first cleared. Then the parameters for the numerical solution of Eq () are defined numerically as in Table 2. Note that a grid in \(x\) is set up as \(x=0,\,0.05,\,0.10,\,\dots,\,1\ .\)
clear all; clc D=1; mu=1; A=1; B=1; nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11; nt=50; nout=10; h=0.001; t=0;
- IC () is then implemented in a for loop (note that t=0).
for i=1:nx; w(i)=A*exp(1/D*(mu^2*t-mu*x(i)))+B; end
- Three nested for loops step Eq () through \(x\) and \(t\)
for i1=1:nout
for i2=1:nt
for i=1:nx
Specifically, the outer loop (in it) gives the nout = 10 outputs in Table 1. Within each output interval (0.05), 50 steps in \(t\) are taken by the loop in i2 (in this way, the 500 steps in \(t\) produce outputs at only 10 points as displayed in Table 1). Within the loop in i, \(x\) spans the interval \(0\leq x\leq 1\ .\)
- BC () is programmed (at i=1) as
if(i==1)w(1)=A*exp(1/D*(mu^2*t-mu*x(1)))+B; wt(1)=0;
Also, since \(w(x=0,t)\) is specified by BC (), the \(t\) derivative of this boundary value is set to zero so that the integration in \(t\) at \(x=0\) does not move \(w(x=0,t)\) from its prescribed value, i.e., wt(1)=0;).
- BC () is used to calculate the fictitious point wf which is then used in the finite difference calculation of Eq ().
elseif(i==nx)wf=w(nx-1)+(2*dx)*A*(-mu/D)*exp(1/D*(mu^2*t-mu*x(nx)));
wt(nx)=D*(wf-2*w(nx)+w(nx-1))/dx^2;
- For the interior points in \(x\) (\(x \neq 0, 1\)), Eq () is programmed.
else wt(i)=D*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
end
The second end terminates the inner loop in i, so that all of the values of \(x\) are handled.
- Euler's method, Eq (), is then used to advance the solution through nt=50 steps of length \(h=0.001\) through the intermediate for loop in i2.
w=w+wt*h; t=t+h; end
Note the use of the matrix facility of MATLAB (w and wt are vectors of length nx = 21). The end completes the intermediate loop in i2.
- The numerical solution is displayed to produce the output in Table 1, including the analytical solution of Eq ()
fprintf('t = %7.4f w(x=xl/2,t) = %7.4f wa(x=xl/2,t) = %7.4f\n', ...
t, w(nx2), A*exp(1/D*(mu^2*t-mu*x(nx2)))+B);
end
The end terminates the outer loop in i1.
Appendix 2
The MATLAB program that produced the results described in Section 4.2 for hyperbolic PDE () (and Eq ()) follows.
clear all; clc
c=1;
nx=201; xl=50; dx=0.5; x=[-xl:dx:xl];
nt=4000; nout=3; h=0.0025; t=0;
for i=1:nx; w(i)=exp(-0.05*x(i)^2); wt(i)=0; ...
wplot(1,i)=w(i); waplot(1,i)=exp(-0.05*x(i)^2); end
for i1=1:nout
for i2=1:nt
for i=1:nx
if(i==1) w( 1)=0; wt( 1)=0; wtt( 1)=0;
elseif(i==nx)w(nx)=0; wt(nx)=0; wtt(nx)=0;
else wtt(i)=c^2*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
end
w=w+wt*h; wt=wt+wtt*h; t=t+h;
end
for i=1:nx
wplot (i1+1,i)=w(i);
waplot(i1+1,i)=0.5*(exp(-0.05*(x(i)-t)^2)+exp(-0.05*(x(i)+t)^2));
fprintf('\n %6.2f %6.2f,%8.4f %8.4f',t, x(i),wplot(i1+1,i),waplot(i1+1,i));
end
end
plot(x,wplot,'-',x,waplot,'o'); xlabel('x'); ylabel('w(x,t)');
title('w(x,t), t=0,10,20,30, solid - num, o - anal')
Listing 2: MATLAB program for the numerical solution of Eqs (), () and ()
We can note the following points about this program:
- Any previous files are first cleared. Then the parameters for the numerical solution of Eq () (or Eq ()) are defined numerically as in Table 3. Note that a grid in \(x\) is set up as \(x=-50,-49.5,\dots,50\ ,\) a total of 201 points. The choice of the upper and lower limits for \(x\ ,\) was made so that domain in \(x\) is essentially infinite (it accommodates BCs () as discussed subsequently).
clear all; clc c=1; nx=201; xl=50; dx=0.5; x=[-xl:dx:xl]; nt=4000; nout=3; h=0.0025; t=0;
- ICs () are then implemented in a for loop (note that t=0).
for i=1:nx; w(i)=exp(-0.05*x(i)^2); wt(i)=0; ...
wplot(1,i)=w(i); waplot(1,i)=exp(-0.05*x(i)^2); end
Both \(w(x,t)\) and \(\dfrac{\partial w(x,t)}{\partial t}\) are programmed (as w and wt) which is required since Eq () is second order in \(t\) (as explained subsequently). Also, the initial value of \(w(x,t=0)\ ,\) both numerical and analytical, is stored for subsequent plotting (so that the plot will include IC ()).
- Three nested for loops step Eqs () and () through \(x\) and \(t\)
for i1=1:nout
for i2=1:nt
for i=1:nx
Specifically, the outer loop (in it) gives the nout = 3 outputs in Table 3. Within each output interval (1), 4000 steps in \(t\) are taken by the loop in i2 (in this way, the 4000 steps in \(t\) produce outputs at only 3 points as displayed in Figure 1). Within the loop in i, \(x\) spans the interval \(-50\leq x\leq 50\ .\)
- BCs () are programmed (at i=1 and i=nx corresponding to \(x=\pm\infty\)) as
if(i==1) w( 1)=0; wt( 1)=0; wtt( 1)=0;
elseif(i==nx)w(nx)=0; wt(nx)=0; wtt(nx)=0;
Also, since \(w\) is specified by BCs (), the first and second derivatives in \(t\) of these boundary values are set to zero so that the integration in \(t\) does not move \(w\) from its zero values.
- For the interior points in \(x\) (\(x \neq -50, 50\)), Eq () (discretized in \(x\)) is programmed.
else wtt(i)=c^2*(w(i+1)-2*w(i)+w(i-1))/dx^2; end
end
The second end terminates the inner loop in i, so that all of the values of \(x\) are handled.
- Euler's method (Eq ()) is then used to advance the solution through nt=4000 steps of length \(h=0.0025\) through the intermediate for loop in i2 according to Eqs () and ().
w=w+wt*h; wt=wt+wtt*h; t=t+h; end
Note the use of the matrix facility of MATLAB (w, wt and wtt are vectors of length nx = 201). The end completes the intermediate loop in i2.
- The numerical and analytical solutions are then stored at \(t=10,20,30\) (three passes through the outer loop in i1) followed by plotting (of Figure 1).
for i=1:nx
wplot (i1+1,i)=w(i);
waplot(i1+1,i)=0.5*(exp(-0.05*(x(i)-t)^2)+exp(-0.05*(x(i)+t)^2));
fprintf('\n %6.2f %6.2f,%8.4f %8.4f',t, x(i),wplot(i1+1,i),waplot(i1+1,i));
end
end
plot(x,wplot,'-',x,waplot,'o'); xlabel('x'); ylabel('w(x,t)');
title('w(x,t), t=0,10,20,30, solid - num, o - anal')
The final end terminates the outer loop in i1. Also, as a word of caution, the fprintf statement displays the complete numerical solution at \(t=10,20,30\) so that the numerical solution can be examined in detail (the peak values reported in Table 4 were taken from this output); if the program in Listing 2 is executed, the fprintf statement could be converted to a comment to reduce the numerical output.
Appendix 3
The MATLAB program that produced the results described in Section 3 for elliptic PDE () (and Eq ()) follows.
clear all; clc
nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11; mu=pi/xl;
ny=21; yl=1; dy=0.05; y=[0:0.05:yl]; ny2=11;
nt=200; nout=10; h=0.0005; t=0;
w=ones(nx,ny);
for i1=1:nout
for i2=1:nt
for i=1:nx
for j=1:ny
if(i== 1) w( 1,j)=0; wt( 1,j)=0;
elseif(i==nx)w(nx,j)=0; wt(nx,j)=0;
elseif(j== 1)wt( i,1)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
+2*(w(i,2)-w(i,1))/dy^2;
elseif(j==ny)wf=w(i,ny-1)+2*dy*mu*sin(mu*x(i))*sinh(mu*y(ny));
wt(i,ny)=(w(i+1,ny)-2*w(i,ny)+w(i-1,ny))/dx^2 ...
+(wf-2*w(i,ny)+w(i,ny-1))/dy^2;
else wt(i,j)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
+(w(i,j+1)-2*w(i,j)+w(i,j-1))/dy^2;
end
end
end
w=w+wt*h;
t=t+h;
end
fprintf('t = %5.2f w(xl/2,yl/2,t) = %7.4f \n',t,w(nx2,ny2))
end
fprintf('\n wa(xl/2,yl/2) = %7.4f\n',sin(mu*x(nx2))*cosh(mu*y(ny2)));
<samll>Listing 3: MATLAB program for the numerical solution of Eqs (), (), (), and ()
We can note the following points about this program:
- Any previous files are first cleared. Then the parameters for the numerical solution of Eq () are defined numerically as in Table 6. Note grids in \(x\) and \(y\) are set up as \(x=0,0.05,0.10,\dots,1\ ,\)\text{ }\(y=0,0.05,0.10,\dots,1\ .\)
clear all; clc nx=21; xl=1; dx=0.05; x=[0:0.05:xl]; nx2=11; mu=pi/xl; ny=21; yl=1; dy=0.05; y=[0:0.05:yl]; ny2=11; nt=200; nout=10; h=0.0005; t=0;
- IC () is then implemented as
w=ones(nx,ny);
\(k=1\) in Eq () is programmed with the MATLAB utility ones over the entire grid in \(x\) and \(y\ ,\) a total of nx \(\times\) ny points.
- Four nested for loops step Eq () through \(x\ ,\) \(y\) and \(t\)
for i1=1:nout
for i2=1:nt
for i=1:nx
for j=1:ny
Specifically, the outer loop (in it) gives the nout = 10 outputs in Table 5. Within each output interval (0.0005), 200 steps in \(t\) are taken by the loop in i2 (in this way, the 2000 steps in \(t\) produce outputs at only 10 points as displayed in Table 5). Within the loop in i, \(x\) spans the interval \(0\leq x\leq 1\) and within the loop in j, \(y\) spans the interval \(0\leq y\leq 1\ .\)
- BC () is programmed (at i=1 for j=1 to j=ny) as
if(i== 1) w( 1,j)=0; wt( 1,j)=0;
Also, since \(w(x=0,y,t)\) is specified by BC (), the \(t\) derivative of this boundary value is set to zero so that the integration in \(t\) at \(x=0\) does not move \(w(x=0,y,t)\) from its prescribed value, i.e., wt(1,j)=0;).
- BC () is programmed (at i=nx for j=1 to j=ny) as
elseif(i==nx)w(nx,j)=0; wt(nx,j)=0;
- Equation () (with subsequent application of Euler's method) is programmed using BC () (at i=1 to i=nx for j=1) as
elseif(j== 1)wt( i,1)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
+2*(w(i,2)-w(i,1))/dy^2;
Note that we have used a finite difference approximation of Eq (), w(i,0) = w(i,2), to eliminate the fictitious value w(i,0) in Eq ().
- Equation () (with subsequent application of Eulers'method) is programmed using BC () (at i=1 to i=nx for j=ny) as
elseif(j==ny)wf=w(i,ny-1)+2*dy*mu*sin(mu*x(i))*sinh(mu*y(ny));
wt(i,ny)=(w(i+1,ny)-2*w(i,ny)+w(i-1,ny))/dx^2 ...
+(wf-2*w(i,ny)+w(i,ny-1))/dy^2;
Note that we have used a finite difference approximation of Eq () to eliminate the fictitious value w(i,ny+1)=wf in Eq ().
- For the interior points in \(x\) (\(x \neq 0, 1\)) and \(y\) (\(x \neq 0, 1\)), Eq ()
(with subsequent application of Euler's method) is programmed.
else wt(i,j)=(w(i+1,j)-2*w(i,j)+w(i-1,j))/dx^2 ...
+(w(i,j+1)-2*w(i,j)+w(i,j-1))/dy^2;
end
end
end
The second and third end statements terminate the inner loops in i and j, so that all of the values of \(x\) and \(y\) are handled.
- Euler's method (Eq ()) is then used to advance the solution through nt=200 steps of length \(h=0.0005\) through the intermediate for loop in i2.
w=w+wt*h;
t=t+h;
end
Note the use of the matrix facility of MATLAB (w and wt are arrays of dimensions nx = 21, ny = 21). The end completes the intermediate loop in i2.
\item The numerical solution is displayed to produce the output in Table 5.
fprintf('t = %5.2f w(xl/2,yl/2,t) = %7.4f \n',t,w(nx2,ny2))
end
The end terminates the outer loop in i1.
- The analytical solution, Eq (), is then printed at the end (as reflected
in Table 5).
fprintf('\n wa(xl/2,yl/2) = %7.4f\n',sin(mu*x(nx2))*cosh(mu*y(ny2)));
The programs in Listings 1, 2 and 3 have an essential feature that should be pointed out - only the operations of storing and retrieving numbers, arithmetic and looping are used (which are done naturally and very efficiently by computers). In other words, the integration of the three PDEs was not done by the usual analytical mathematics, but rather, was done essentially with arithmetic. This is the essence of the numerical method, that is, to replace mathematics of differentiation and integration with something much simpler, but performed naturally by computers with extraordinary speed and precision. This simplification was accomplished by replacing the partial derivatives in the PDEs with algebraic approximations, in this case, finite differences.
A central question then is how accurate the numerical solution is if approximations are used to compute it. We have observed the accuracy can be improved with more calculations (e.g., by reducing \(\Delta x\ ,\) \(\Delta y\) and \(h\)). In other words there is a tradeoff between accuracy and calculations, and we basically assume the computer can do enough calculations to achieve the required accuracy. Historically, the increasing speed of computers has produced solutions to PDE problems with acceptable accuracy and continually increasing complexity. Thus, for the future, we can expect to compute numerical solutions of PDEs with essentially unlimited complexity if sufficient computing power is available.
