[With code and PDF] How can I solve systems of PDEs with three independent variables by using the ode15s solver of the MATLAB?

6 views (last 30 days)
ponqun
ponqun on 27 Jan 2015
Commented: Torsten on 27 Jan 2015
[***Please find the matlab commands and PDF files***]
Now, I am trying to solve systems of PDEs (attached PDF file) by MATLAB. As you can see, this PDEs have three independent variable (z,r,t).
According to a kind suggestion by "Torsten" (<http://www.mathworks.com/matlabcentral/answers/170648)>, I have discretized my equations about z- and r-coordinate. Then, I have developed matlab codes (PDEbook_s.m with pde_s.m, dss004,m and dss044.m) by using ode15s as follows. dss004 and dss044 (by http://www.pdecomp.net/) calculate first and second derivative values, respectively.
In this PDEs, the U and C values should decrease to zero value with increasing time. However, this values closed to a constant value. (Please run the commands in your command window) . Therefore, I think that this command cannot return correct values of the PDE systems now.
Can you find any problem in these codes?
*************PDEbook_s.m ****************************************************************************************
format short
%
% Clear previous files
clear all
clc
%
%Global data
global nr nz dr dz drs dzs...
r z zl r0 call tf ncall...
%
%Model parameters
zl=1;
r0=1;
%
% Grid in axial direcion
nz=10;
dz=zl/(nz-1);
for i=1:nz;
z(i)=(i-1)*dz;
end
dzs=dz^2;
%
%Grid in radial direction
nr=10;
dr=r0/(nr-1);
for j=1:nr
r(j)=(j-1)*dr;
end
drs=dr^2;
%
%Independent variables for ODE integration
tf=5;
delt=1.0E-2;
tout=[0.0:delt:(tf-delt)];
nout=(tf)/delt;
call=0;
%
%intial condition
for i=1:nz
for j=1:nr
C(i,j)=0;
U(i,j)=1;
y0((i-1)*nr+j)=C(i,j);
y0((i-1)*nr+j+nz*nr)=U(i,j);
end
end
%
%ODE integration
reltol=1.0E-4; abstol=1.0E-5;
options=odeset('RelTol',reltol, 'AbsTol',abstol);
[t,y]=ode15s(@pde_s,tout,y0,options);
%
%1D to 2D matrices
for it=1:nout
for i=1:nz
for j=1:nr
C(it,i,j)=y(it,(i-1)*nr+j);
U(it,i,j)=y(it,(i-1)*nr+j+nz*nr);
end
end
end
%
% t-C,U plots (r=r0)
figure(1);
subplot(2,2,1)
plot(tout,C(:,1,nr),'--g',tout,U(:,1,nr),':r'); axis([0 tf 0 1]);
title('Solution(z=0,r=r0)')
xlabel('time')
ylabel('C(green),U(red)')
subplot(2,2,2)
plot(tout,C(:,5,nr),'--g',tout,U(:,5,nr),':r'); axis([0 tf 0 1]);
title('Solution(z=4*dz,r=r0)')
xlabel('time')
ylabel('C(green),U(red)')
subplot(2,2,3)
plot(tout,C(:,8,nr),'--g',tout,U(:,8,nr),':r'); axis([0 tf 0 1]);
title('Solution(z=*7dz,r=r0)')
xlabel('time')
ylabel('C(green),U(red)')
subplot(2,2,4)
plot(tout,C(:,nz,nr),'--g',tout,U(:,nz,nr),':r'); axis([0 tf 0 1]);
title('Solution(z=zl,r=r0)')
xlabel('time')
ylabel('C(green),U(red)')
%
% t-C,U plots (r=1/2r0)
figure(2);
subplot(2,2,1)
plot(tout,C(:,1,5),'--g',tout,U(:,1,5),':r'); axis([0 tf 0 1]);
title('Solution(=0,r=1/2r0)')
xlabel('time')
ylabel('C(green),U(red)')
subplot(2,2,2)
plot(tout,C(:,5,5),'--g',tout,U(:,5,5),':r'); axis([0 tf 0 1]);
title('Solution(=4*dz,r=1/2r0)')
xlabel('time')
ylabel('C(green),U(red)')
subplot(2,2,3)
plot(tout,C(:,8,5),'--g',tout,U(:,8,5),':r'); axis([0 tf 0 1]);
title('Solution(z=*7dz,r=1/2r0)')
xlabel('time')
ylabel('C(green),U(red)')
subplot(2,2,4)
plot(tout,C(:,nz,5),'--g',tout,U(:,nz,5),':r'); axis([0 tf 0 1]);
title('Solution(z=zl,r=1/2r0)')
xlabel('time')
ylabel('C(green),U(red)')
*****End of PDEbook_s****************************************************************************************
***********pde_s.m ****************************************************************************************
function yt=pde_s(t,y)
%
%Global data
global nr nz dr dz drs dzs...
r z zl r0 call tf ncall...
%
% 1D to 2D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
C(i,j)=y(ij);
U(i,j)=y(ij+nr*nz);
end
end
%
% Step through the grid points in r and z
for i=1:nz
C1d=C(i,:); %1d means 1D value
U1d=U(i,:);
%
% Ur by dss004:Ur means first derivative of U for r
Ur1d=dss004(0.0,r0,nr,U1d); %Ur means first derivative of U for r
Ur(i,:)=Ur1d;
Ur(i,1)= 0.0; %B.C.1 of Eq.(1)
Ur(i,nr)=C(i,j)-U(i,nr);%B.C.2 of Eq.(1)
%
%Urr by dss 044 with the Neumann B.D. (nl=nu=2):Urr means second derivative of U for r
Ur1d( 1)=0.0;%B.C.1 of Eq.(1)
Ur1d(nr)=C1d(j)-U1d(nr);%B.C.2 of Eq.(1)
nl=2;
nu=2;
Urr1d=dss044(0.0,r0,nr,U1d,Ur1d,nl,nu);
Urr(i,:)=Urr1d;
%
for j=1:nr
%
% (1/r)*Ur
if(j~=1)
Ur(i,j)=(1.0/r(j))*Ur(i,j);
end
if(j==1)Urr(i,j)=2.0*Urr(i,j);end
%
% Cz by dss004:Cz means first derivative of C for z
Cz1d=dss004(0.0,zl,nz,C1d);
Cz(:,j)=Cz1d;
Cz(1,j)= C1d(1);%B.C.1 of Eq.(2)
Cz(nz,j)=0.0;%B.C.2 of Eq.(2)
%
% Czz by dss044 with the Neumann B.D. (nl=nu=2):Czz means second derivative of C for z
% Czz
Cz1d(1)=C(1,nr);%B.C.1 of Eq.(2)
Cz1d(nz)=0.0;%B.C.2 of Eq.(2)
nl=2;
nu=2;
Czz1d=dss044(0.0,zl,nz,C1d,Cz1d,nl,nu);
Czz(:,j)=Czz1d;
%
% PDEs
Ut(i,j)=(1/(1+U(i,j)))*Ur(i,j)+(1/(1+U(i,j)))*Urr(i,j); %Eq.(1)
Ct(i,j)=Czz(i,j)-Cz(i,j)+(U(i,nr)-C(i,j)); % Eq.(2)
end
end
%
% 2D to 1D matrices
for i=1:nz
for j=1:nr
ij=(i-1)*nr+j;
yt(ij)=Ct(i,j);
yt(ij+nr*nz)=Ut(i,j);
end
end
%
% Transpose and count
yt=yt';
ncall=ncall+1;
end
***********End of pde_s.m ****************************************************************************************
***********dss004.m ****************************************************************************************
% File: dss004.m
%
function [ux]=dss004(xl,xu,n,u)
%
% Function dss004 computes the first derivative, u , of a
% x
% variable u over the spatial domain xl le x le xu from classical
% five-point, fourth-order finite difference approximations
%
% Argument list
%
% xl Lower boundary value of x (input)
%
% xu Upper boundary value of x (input)
%
% n Number of grid points in the x domain including the
% boundary points (input)
%
% u One-dimensional array containing the values of u at
% the n grid point points for which the derivative is
% to be computed (input)
%
% ux One-dimensional array containing the numerical
% values of the derivatives of u at the n grid points
% (output)
%
% The mathematical details of the following Taylor series (or
% polynomials) are given in routine dss002.
%
% Five-point formulas
%
% (1) Left end, point i = 1
%
% 2 3 4
% a(u2 = u1 + u1 ( dx) + u1 ( dx) + u1 ( dx) + u1 ( dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 ( dx) + u1 ( dx) + u1 ( dx) + ...)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% b(u3 = u1 + u1 (2dx) + u1 (2dx) + u1 (2dx) + u1 (2dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (2dx) + u1 (2dx) + u1 (2dx) + ...)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% c(u4 = u1 + u1 (3dx) + u1 (3dx) + u1 (3dx) + u1 (3dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (3dx) + u1 (3dx) + u1 (3dx) + ...)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% d(u5 = u1 + u1 (4dx) + u1 (4dx) + u1 (4dx) + u1 (4dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u1 (4dx) + u1 (4dx) + u1 (4dx) + ...)
% 5x 5f 6x 6f 7x 7f
%
% Constants a, b, c and d are selected so that the coefficients
% of the u1 terms sum to one and the coefficients of the u1 ,
% x 2x
% u1 and u1 terms sum to zero
% 3x 4x
%
% a + 2b + 3c + 4d = 1
%
% a + 4b + 9c + 16d = 0
%
% a + 8b + 27c + 64d = 0
%
% a + 16b + 81c + 256d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% u1 = (1/12dx)(-25u1 + 48u2 - 36u3 + 16u4 - 3u5) + O(dx ) (1)
% x
%
% (2) Interior point, i = 2
%
% 2 3 4
% a(u1 = u2 + u2 (-dx) + u2 (-dx) + u2 (-dx) + u2 (-dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (-dx) + u2 (-dx) + u2 (-dx) + ...)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% b(u3 = u2 + u2 ( dx) + u2 ( dx) + u2 ( dx) + u2 ( dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 ( dx) + u2 ( dx) + u2 ( dx) + ...)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% c(u4 = u2 + u2 (2dx) + u2 (2dx) + u2 (2dx) + u2 (2dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (2dx) + u2 (2dx) + u2 (2dx) + ...)
% 5x 5f 6x 6f 7x 7f
%
% 2 3 4
% d(u5 = u2 + u2 (3dx) + u2 (3dx) + u2 (3dx) + u2 (3dx)
% x 1f 2x 2f 3x 3f 4x 4f
%
% 5 6 7
% + u2 (3dx) + u2 (3dx) + u2 (3dx) + ...)
% 5x 5f 6x 6f 7x 7f
%
% -a + b + 2c + 3d = 1
%
% a + b + 4c + 9d = 0
%
% -a + b + 8c + 27d = 0
%
% a + b + 16c + 81d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% u2 = (1/12dx)(-3u1 - 10u2 + 18u3 - 6u4 + u5) + O(dx ) (2)
% x
%
% (3) Interior point i, i ne 2, n-1
%
% 2 3
% a(ui-2 = ui + ui (-2dx) + ui (-2dx) + ui (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui (-2dx) + ui (-2dx) + ui (-2dx) + ...)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(ui-1 = ui + ui ( -dx) + ui ( -dx) + ui ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( -dx) + ui ( -dx) + ui ( -dx) + ...)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(ui+1 = ui + ui ( dx) + ui ( dx) + ui ( dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( dx) + ui ( dx) + ui ( dx) + ...)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(ui+2 = ui + ui ( 2dx) + ui ( 2dx) + ui ( 2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + ui ( 2dx) + ui ( 2dx) + ui ( 2dx) + ...)
% 4x 4f 5x 5f 6x 6f
%
% -2a - b + c + 2d = 1
%
% 4a + b + c + 4d = 0
%
% -8a - b + c + 8d = 0
%
% 16a + b + c + 16d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% ui = (1/12dx)(ui-2 - 8ui-1 + 0ui + 8ui+1 - ui+2) + O(dx ) (3)
% x
%
% (4) Interior point, i = n-1
%
% 2 3
% a(un-4 = un-1 + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 (-3dx) + un-1 (-3dx) + un-1 (-3dx) + ...
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(un-3 = un-1 + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 (-2dx) + un-1 (-2dx) + un-1 (-2dx) + ...
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(un-2 = un-1 + un-1 ( -dx) + un-1 (- -x) + un-1 ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 ( -dx) + un-1 ( -dx) + un-1 ( -dx) + ...
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(un = un-1 + un-1 ( dx) + un-1 ( dx) + un-1 ( dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un-1 ( dx) + un-1 ( dx) + un-1 ( dx) + ...
% 4x 4f 5x 5f 6x 6f
%
% -3a - 2b - c + d = 1
%
% 9a + 4b + c + d = 0
%
% -27a - 8b - c + d = 0
%
% 81a + 16b + c + d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% un-1 = (1/12dx)(-un-4 + 6un-3 - 18un-2 + 10un-1 + 3un) + O(dx )
% x
% (4)
%
% (5) Right end, point i = n
%
% 2 3
% a(un-4 = un + un (-4dx) + un (-4dx) + un (-4dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-4dx) + un (-4dx) + un (-4dx) + ...)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% b(un-3 = un + un (-3dx) + un (-3dx) + un (-3dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-3dx) + un (-3dx) + un (-3dx) + ...)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% c(un-2 = un + un (-2dx) + un (-2dx) + un (-2dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un (-2dx) + un (-2dx) + un (-2dx) + ...)
% 4x 4f 5x 5f 6x 6f
%
% 2 3
% d(un-1 = un + un ( -dx) + un ( -dx) + un ( -dx)
% x 1f 2x 2f 3x 3f
%
% 4 5 6
% + un ( -dx) + un ( -dx) + un ( -dx) + ...)
% 4x 4f 5x 5f 6x 6f
%
% -4a - 3b - 2c - d = 1
%
% 16a + 9b + 4c + d = 0
%
% -64a - 27b - 8c - d = 0
%
% 256a + 81b + 16c + d = 0
%
% Simultaneous solution for a, b, c and d followed by the solu-
% tion of the preceding Taylor series, truncated after the u
% 4x
% terms, for u1 gives the following five-point approximation
% x
% 4
% un = (1/12dx)(3un-4 - 16un-3 + 36un-2 - 48un-1 + 25un) + O(dx )
% x
% (5)
%
% The weighting coefficients for equations (1) to (5) can be
% summarized as
%
% -25 48 -36 16 -3
%
% -3 -10 18 -6 1
%
% 1/12 1 -8 0 8 -1
%
% -1 6 -18 10 3
%
% 3 -16 36 -48 25
%
% which are the coefficients reported by Bickley for n = 4, m =
% 1, p = 0, 1, 2, 3, 4 (Bickley, W. G., Formulae for Numerical
% Differentiation, Math. Gaz., vol. 25, 1941. Note - the Bickley
% coefficients have been divided by a common factor of two).
%
% Equations (1) to (5) can now be programmed to generate the
% derivative u (x) of function u(x).
% x
%
% Compute the spatial increment
dx=(xu-xl)/(n-1);
r4fdx=1./(12.*dx);
nm2=n-2;
%
% Equation (1) (note - the rhs of equations (1), (2), (3), (4)
% and (5) have been formatted so that the numerical weighting
% coefficients can be more easily associated with the Bickley
% matrix above)
ux( 1)=r4fdx*...
( -25.*u( 1) +48.*u( 2) -36.*u( 3) +16.*u( 4) -3.*u( 5));
%
% Equation (2)
ux( 2)=r4fdx*...
( -3.*u( 1) -10.*u( 2) +18.*u( 3) -6.*u( 4) +1.*u( 5));
%
% Equation (3)
for i=3:nm2
ux( i)=r4fdx*...
( +1.*u(i-2) -8.*u(i-1) +0.*u( i) +8.*u(i+1) -1.*u(i+2));
end
%
% Equation (4)
ux(n-1)=r4fdx*...
( -1.*u(n-4) +6.*u(n-3) -18.*u(n-2) +10.*u(n-1) +3.*u( n));
%
% Equation (5)
ux( n)=r4fdx*...
( 3.*u(n-4) -16.*u(n-3) +36.*u(n-2) -48.*u(n-1) +25.*u( n));
***********End of dss004.m **************************************************************************************
***********dss044.m ****************************************************************************************
% File: ds044.m
%
function [uxx]=dss044(xl,xu,n,u,ux,nl,nu)
%
% Function dss044 computes a fourth-order approximation of a
% second-order derivative, with or without the normal derivative
% at the boundary.
%
% Argument list
%
% xl Left value of the spatial independent variable (input)
%
% xu Right value of the spatial independent variable (input)
%
% n Number of spatial grid points, including the end
% points (input)
%
% u One-dimensional array of the dependent variable to be
% differentiated (input)
%
% ux One-dimensional array of the first derivative of u.
% The end values of ux, ux(1) and ux(n), are used in
% Neumann boundary conditions at x = xl and x = xu,
% depending on the arguments nl and nu (see the de-
% scription of nl and nu below)
%
% uxx One-dimensional array of the second derivative of u
% (output)
%
% nl Integer index for the type of boundary condition at
% x = xl (input). The allowable values are
%
% 1 - Dirichlet boundary condition at x = xl
% (ux(1) is not used)
%
% 2 - Neumann boundary condition at x = xl
% (ux(1) is used)
%
% nu Integer index for the type of boundary condition at
% x = xu (input). The allowable values are
%
% 1 - Dirichlet boundary condition at x = xu
% (ux(n) is not used)
%
% 2 - Neumann boundary condition at x = xu
% (ux(n) is used)
%
% The following derivation was completed by W. E. Schiesser, Depts
% of CHE and Math, Lehigh University, Bethlehem, PA 18015, USA, on
% December 15, 1986. Additional details are given in function
% dss042.
%
% ******************************************************************
%
% (1) uxx at the interior points 3, 4,..., n-2
%
% To develop a set of fourth-order correct differentiation formulas
% for the second derivative uxx, we consider first the interior
% grid points at which a symmetric formula can be used.
%
% If we consider a formula of the form
%
% a*u(i-2) + b*u(i-1) + e*u(i) + c*u(i+1) + d*u(i+2)
%
% Taylor series expansions of u(i-2), u(i-1), u(i+1) and u(i+2)
% can be substituted into this formula. We then consider the
% linear albegraic equations relating a, b, c and d which will
% retain certain terms, i.e., uxx, and drop others, e.g., uxxx,
% uxxxx and uxxxxx.
%
% Thus, for grid points 3, 4,..., n-2
%
% To retain uxx
%
% 4*a + b + c + 4*d = 2 (1)
%
% To drop uxxx
%
% -8*a - b + c + 8*d = 0 (2)
%
% To drop uxxxx
%
% 16*a + b + c + 16*d = 0 (3)
%
% To drop uxxxxx
%
% -32*a - b + c + 32*d = 0 (4)
%
% Equations (1) to (4) can be solved for a, b, c and d. If equa-
% tion (1) is added to equation (2)
%
% -4*a + 2*c + 12*d = 2 (5)
%
% If equation (1) is subtracted from equation (3)
%
% 12*a + 12*d = -2 (6)
%
% If equation (1) is added to equation (4)
%
% -28*a + 2*c + 36*d = 2 (7)
%
% Equations (5) to (7) can be solved for a, c and d. If equation
% (5) is subtracted from equation (7), and the result combined
% with equation (6)
%
% 12*a + 12*d = -2 (6)
%
% -24*a + 24*d = 0 (8)
%
% Equations (6) and (8) can be solved for a and d. From (8), a
% = d. From equation (6), a = -1/12 and d = -1/12. Then, from
% equation (5), c = 4/3, and from equation (1), b = 4/3.
%
% The final differentiation formula is then obtained as
%
% (-1/12)*u(i-2) + (4/3)*u(i-1) +
%
% (4/3)*u(i+1) + (-1/12)*u(i+2)
%
% (-1/12 + 4/3 - 1/12 + 4/3)*u(i) + uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/(12*dx**2))*(-1*u(i-2) + 16*u(i-1)
%
% - 30*u(i) + 16*u(i+1) - 1*u(i+2) (9)
%
% + O(dx**4)
%
% Note that the ux term drops out, i.e., the basic equation is
%
% -2*a - b + c + 2*d =
%
% -2*(-1/12) - (4/3) + (4/3) + 2*(-1/12) = 0
%
% Equation (9) was obtained by dropping all terms in the underlying
% Taylor series up to and including the fifth derivative, uxxxxx.
% Thus, equation (9) is exact for polynomials up to and including
% fifth order. This can be checked by substituting the functions
% 1, x, x**2, x**3, x**4 and x**5 in equation (9) and computing the
% corresponding derivatives for comparison with the known second
% derivatives. This is done for 1 merely by summing the weighting
% coefficients in equation (9), which should sum to zero, i.e.,
% -1 + 16 - 30 + 16 -1 = 0.
%
% For the remaining functions, the algebra is rather involved, but
% these functions can be checked numerically, i.e., numerical values
% of x**2, x**3, x**4 and x**5 can be substituted in equation (9)
% and the computed derivatives can be compared with the know numeri-
% cal second derivatives. This is not a proof of correctness of
% equation (9), but would likely detect any errors in equation (9).
%
% ******************************************************************
%
% (2) uxx at the interior points i = 2 and n-1
%
% For grid point 2, we consider a formula of the form
%
% a*u(i-1) + f*u(i) + b*u(i+1) + c*u(i+2) + d*u(i+3) + e*u(i+4)
%
% Taylor series expansions of u(i-1), u(i+1), u(i+2), u(i+3) and
% u(i+4) when substituted into this formula give linear algebraic
% equations relating a, b, c, d and e.
%
% To drop ux
%
% -a + b + 2*c + 3*d + 4*e = 0 (10)
%
% To retain uxx
%
% a + b + 4*c + 9*d + 16*e = 2 (11)
%
% To drop uxxx
%
% -a + b + 8*c + 27*d + 64*e = 0 (12)
%
% To drop uxxxx
%
% a + b + 16*c + 81*d + 256*e = 0 (13)
%
% To drop uxxxxx
%
% -a + b + 32*c + 243*d + 1024*e = 0 (14)
%
% Equations (11), (12), (13) and (14) can be solved for a, b, c,
% d and e. If equation (10) is added to equation (11)
%
% 2*b + 6*c + 12*d +20*e = 2 (15)
%
% If equation (10) is subtracted from equation (12)
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (10) is added to equation (13)
%
% 2*b + 18*c + 84*d + 260*e = 0 (17)
%
% If equation (10) is subtracted from equation (14)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (15), (16), (17) and (18) can be solved for b, c, d
% and e.
%
% 6*c + 24*d + 60*e = 0 (16)
%
% If equation (15) is subtracted from equation (17)
%
% 12*c + 72*d + 240*e = -2 (19)
%
% 30*c + 240*d + 1020*e = 0 (18)
%
% Equations (16), (18) and (19) can be solved for c, d and e. If
% two times equation (16) is subtracted from equation (19),
%
% 24*d + 120*e = -2 (20)
%
% If five times equation (16) is subtracted from equation (18),
%
% 120*d + 720*e = 0 (21)
%
% Equations (20) and (21) can be solved for d and e. From (21),
% e = (-1/6)*d. Substitution in equation (20) gives d = -1/2.
% thus, e = 1/12. From equation (16), c = 7/6. From equation
% (15), b = -1/3. From equation (10), a = 5/6.
%
% The final differentiation formula is then obtained as
%
% (5/6)*u(i-1) + (-1/3)*u(i+1) + (7/6)*u(i+2) + (-1/2)*u(i+3)
%
% + (1/12)*u(i+4) = (5/6 - 1/3 + 7/6 - 1/2 + 1/12)*u(i)
%
% + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(10*u(i-1) - 15*u(i) - 4*u(i+1)
% (22)
% + 14*u(i+2) - 6*u(i+3) + 1*u(i+4)) + O(dx**4)
%
% Equation (22) will be applied at i = 2 and n-1. thus
%
% uxx(2) = (1/12*dx**2)*(10*u(1) - 15*u(2) - 4*u(3)
% (23)
% + 14*u(4) - 6*u(5) + 1*u(6)) + O(dx**4)
%
% uxx(n-1) = (1/12*dx**2)*(10*u(n) - 15*u(n-1) - 4*u(n-2)
% (24)
% + 14*u(n-3) - 6*u(n-4) + 1*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% (3) uxx at the boundary points 1 and n
%
% Finally, for grid point 1, an approximation with a Neumann bound-
% ary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*ux(i) + f*u(i)
%
% Will be used. the corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + e = 0 (25)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d = 2 (26)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d = 0 (27)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d = 0 (28)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d = 0 (29)
%
% Equations (25) to (29) can be solved for a, b, c, d and e. If
%
% Equation (26) is subtracted from equations (27), (28) and (29),
%
% 4*b + 18*c + 48*d = -2 (30)
%
% 12*b + 72*c + 240*d = -2 (31)
%
% 28*b + 234*c + 1008*d = -2 (32)
%
% Equations (30), (31) and (32) can be solved for b, c and d
%
% 18*c + 96*d = 4 (33)
%
% 108*c + 672*d = 12 (34)
%
% Equations (3) and (34) can be solved for c and d, c = 8/9, d =
% -1/8.
%
% From equation (30), b = -3. From equation (26), a = 8. From
% equation (25), e = -25/6.
%
% The final differentiation formula is then obtained as
%
% 8*u(i+1) - 3*u(i+2) + (8/9)*u(i+3) - (1/8)*u(i+4)
%
% - (25/6)*ux(i)*dx
%
% = (8 - 3 + (8/9) - (1/8))*u(i) + uxx*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*((-415/6)*u(i) + 96*u(i+1) - 36*u(i+2)
% (35)
% + (32/3)*u(i+3) - (3/2)*u(i+4) - 50*ux(i)*dx) + O(dx**4)
%
% Equation (35) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*((-415/6)*u(1) + 96*u(2) - 36*u(3)
% (36)
% + (32/3)*u(4) - (3/2)*u(5) - 50*ux(1)*dx) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*((-415/6)*u(n) + 96*u(n-1) - 36*u(n-2)
% (37)
% + (32/3)*u(n-3) - (3/2)*u(n-4) + 50*ux(n)*dx) + O(dx**4)
%
% Alternatively, for grid point 1, an approximation with a Dirichlet
% boundary condition of the form
%
% a*u(i+1) + b*u(i+2) + c*u(i+3) + d*u(i+4) + e*u(i+5) + f*u(i)
%
% can be used. The corresponding algebraic equations are
%
% To drop ux
%
% a + 2*b + 3*c + 4*d + 5*e = 0 (38)
%
% To retain uxx
%
% a + 4*b + 9*c + 16*d + 25*e = 2 (39)
%
% To drop uxxx
%
% a + 8*b + 27*c + 64*d + 125*e = 0 (40)
%
% To drop uxxxx
%
% a + 16*b + 81*c + 256*d + 625*e = 0 (41)
%
% To drop uxxxxx
%
% a + 32*b + 243*c + 1024*d + 3125*e = 0 (42)
%
% Equations (38), (39), (40), (41) amd (42) can be solved for a,
% b, c, d and e.
%
% 2*b + 6*c + 12*d + 20*e = 2 (43)
%
% 6*b + 24*c + 60*d + 120*e = 0 (44)
%
% 14*b + 78*c + 252*d + 620*e = 0 (45)
%
% 30*b + 240*c + 1020*d + 3120*e = 0 (46)
%
% Equations (43), (44), (45) and (46) can be solved for b, c, d
% and e
%
% 6*c + 24*d + 60*e = -6 (47)
%
% 36*c + 168*d + 480*e = -14 (48)
%
% 150*c + 840*d + 2820*e = -30 (49)
%
% Equations (47), (48) and (49) can be solved for c, d and e
%
% 24*d + 120*e = 22 (50)
%
% 240*d + 1320*e = 120 (51)
%
% From equations (50) and (51), d = 61/12, e = -5/6. From equation
% (47), c = -13. From equation (43), b = 107/6. From equation (38),
% a = -77/6.
%
% The final differentiation formula is then obtained as
%
% (-77/6)*u(i+1) + (107/6)*u(i+2) - 13*u(i+3) + (61/12)*u(i+4)
%
% - (5/6)*u(i+5) = (-77/6 + 107/6 - 13 + 61/12 - 5/6)*u(i) +
%
% uxx(i)*(dx**2) + O(dx**6)
%
% or
%
% uxx(i) = (1/12*dx**2)*(45*u(i) - 154*u(i+1) + 214*u(i+2)
% (52)
% - 156*u(i+3) + 61*u(i+4) - 10*u(i+5)) + O(dx**4)
%
% Equation (52) will be applied at i = 1 and i = n
%
% uxx(1) = (1/12*dx**2)*(45*u(1) - 154*u(2) + 214*u(3)
% (53)
% - 156*u(4) + 61*u(5) - 10*u(6)) + O(dx**4)
%
% uxx(n) = (1/12*dx**2)*(45*u(n) - 154*u(n-1) + 214*u(n-2)
% (54)
% -156*u(n-3) + 61*u(n-4) - 10*u(n-5)) + O(dx**4)
%
% ******************************************************************
%
% Grid spacing
dx=(xu-xl)/(n-1);
%
% 1/(12*dx**2) for subsequent use
r12dxs=1./(12.0*dx^2);
%
% uxx at the left boundary
%
% Without ux (equation (53))
if nl==1
uxx(1)=r12dxs*...
( 45.0*u(1)...
-154.0*u(2)...
+214.0*u(3)...
-156.0*u(4)...
+61.0*u(5)...
-10.0*u(6));
%
% With ux (equation (36))
elseif nl==2
uxx(1)=r12dxs*...
(-415.0/6.0*u(1)...
+96.0*u(2)...
-36.0*u(3)...
+32.0/3.0*u(4)...
-3.0/2.0*u(5)...
-50.0*ux(1)*dx);
end
%
% uxx at the right boundary
%
% Without ux (equation (54))
if nu==1
uxx(n)=r12dxs*...
( 45.0*u(n )...
-154.0*u(n-1)...
+214.0*u(n-2)...
-156.0*u(n-3)...
+61.0*u(n-4)...
-10.0*u(n-5));
%
% With ux (equation (37))
elseif nu==2
uxx(n)=r12dxs*...
(-415.0/6.0*u(n )...
+96.0*u(n-1)...
-36.0*u(n-2)...
+32.0/3.0*u(n-3)...
-3.0/2.0*u(n-4)...
+50.0*ux(n )*dx);
end
%
% uxx at the interior grid points
%
% i = 2 (equation (23))
uxx(2)=r12dxs*...
( 10.0*u(1)...
-15.0*u(2)...
-4.0*u(3)...
+14.0*u(4)...
-6.0*u(5)...
+1.0*u(6));
%
% i = n-1 (equation (24))
uxx(n-1)=r12dxs*...
( 10.0*u(n )...
-15.0*u(n-1)...
-4.0*u(n-2)...
+14.0*u(n-3)...
-6.0*u(n-4)...
+1.0*u(n-5));
%
% i = 3, 4,..., n-2 (equation (9))
for i=3:n-2
uxx(i)=r12dxs*...
( -1.0*u(i-2)...
+16.0*u(i-1)...
-30.0*u(i )...
+16.0*u(i+1)...
-1.0*u(i+2));
end
***********End of dss044.m **************************************************************************************
  1 Comment
Torsten
Torsten on 27 Jan 2015
Maybe I'm mistaken, but I guess the number of equations you have to solve is nz*nr for the variable U and nz for the variable C. So your C-array should be one-dimensional, shouldn't it ?
Best wishes
Torsten.

Sign in to comment.

Sign in to answer this question.

Answers (0)

Categories

Find more on General PDEs in Help Center and File Exchange

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!