Boundary conditions with the diffusion operator step
A tutorial of the variable diffusivity diffusion equation based on Script 36 from Trefethen's Book. This uses the LU decomposition. The steps are also timed. This is an extension of the steady_diff_bc script. Please note that this script defines functions at the end, which is only supported by MATLAB 2016b or later.
Contents
clear,close all
Set up grid and 2D Laplacian, boundary points included:
Notice the low number of points.
N = 64;
Use the cheb function to get the Chebyshev points.
[D,x] = cheb(N); y = x;
Create the 2D grid.
[xx,yy] = meshgrid(x,y); [xxx,yyy] = meshgrid(-1:.04:1,-1:.04:1);
Create the differentiation matrices
Make them sparse to save space.
D2 = sparse(D^2);
I = speye(N+1);
% Use the outer product to create the matrices.
Dy=kron(I,D);
Dx=kron(D,I);
Dyy=kron(I,D2);
Dxx=kron(D2,I);
Create variable diffusion constant
r=sqrt((xx-0.5).^2+(yy+0.75).^2); kappa0=1e-2; kappa=kappa0*(1+0.5*sech((r/0.5).^2)); kappav=kappa(:); mykappax=Dx*kappa(:); mykappay=Dy*kappa(:); kappaxv=mykappax(:); kappayv=mykappay(:); kappavmat=sparse(diag(kappav)); kappaxmat=sparse(diag(kappaxv)); kappaymat=sparse(diag(kappayv));
Create the Laplacian operator
Combine the last two parts and build the Laplacian operator.
Llap = kappavmat*(Dxx+Dyy)+kappaxmat*Dx+kappaymat*Dy;
Create the time evolving Laplace operator
dt=0.1; L=speye(size(Llap))-dt*Llap;
Boundary Conditions
Impose boundary conditions by replacing appropriate rows of L:
bxp1 = find(xx==1); % Rows of the big matrix at which x=1. bxn1 = find(xx==-1); % Rows of the big matrix at which x=-1. byp1 = find(yy==1); % Rows of the big matrix at which y=1. byn1 = find(yy==-1); % Rows of the big matrix at which y=-1. rhs = zeros((N+1)^2,1); L(bxp1,:) = zeros(N+1,(N+1)^2); L(bxp1,bxp1) = eye(N+1); rhs(bxp1)=-yy(bxp1); % Dirichlet BC at x=-1 notice how RHS is modified. for ii=1:length(byp1) L(byp1(ii),:) = Dy(byp1(ii),:); rhs(byp1)=2.5; % Set L=Dy at y=1, notice RHS is set to zero above. L(bxn1(ii),:) = Dx(bxn1(ii),:); rhs(bxn1)=sin(pi*(1+yy(bxn1))); % Set L=Dy at x=-1, notice RHS modifed. end L(byn1,:) = zeros(N+1,(N+1)^2); L(byn1,byn1) = eye(N+1); rhs(byn1)=-cos(3*pi*xx(byn1));
Solve the equation
Use spparms to show the logic used by Matlab to solve the linear system with the "\" command.
spparms('spumoni',2)
Just solve the equations with nothing extra.
tic u1 = L\rhs; uu1 = reshape(u1,N+1,N+1); justsolve=toc
sp\: bandwidth = 4160+1+4095.
sp\: is A diagonal? no.
sp\: is band density (0.0291729) > bandden (0.5) to try banded solver? no.
sp\: is A triangular? no.
sp\: is A morally triangular? no.
sp\: is A a candidate for Cholesky (symmetric, real positive diagonal)? no.
sp\: use Unsymmetric MultiFrontal PACKagewith automatic reordering.
UMFPACK V5.4.0 (May 20, 2009), Control:
Matrix entry defined as: double
Int (generic integer) defined as: UF_long
0: print level: 2
1: dense row parameter: 0.2
"dense" rows have > max (16, (0.2)*16*sqrt(n_col) entries)
2: dense column parameter: 0.2
"dense" columns have > max (16, (0.2)*16*sqrt(n_row) entries)
3: pivot tolerance: 0.1
4: block size for dense matrix kernels: 32
5: strategy: 0 (auto)
6: initial allocation ratio: 0.7
7: max iterative refinement steps: 2
12: 2-by-2 pivot tolerance: 0.01
13: Q fixed during numerical factorization: 0 (auto)
14: AMD dense row/col parameter: 10
"dense" rows/columns have > max (16, (10)*sqrt(n)) entries
Only used if the AMD ordering is used.
15: diagonal pivot tolerance: 0.001
Only used if diagonal pivoting is attempted.
16: scaling: 1 (divide each row by sum of abs. values in each row)
17: frontal matrix allocation ratio: 0.5
18: drop tolerance: 0
19: AMD and COLAMD aggressive absorption: 1 (yes)
The following options can only be changed at compile-time:
8: BLAS library used: Fortran BLAS. size of BLAS integer: 8
9: compiled for MATLAB
10: CPU timer is POSIX times ( ) routine.
11: compiled for normal operation (debugging disabled)
computer/operating system: Linux
size of int: 4 UF_long: 8 Int: 8 pointer: 8 double: 8 Entry: 8 (in bytes)
sp\: UMFPACK's factorization was successful.
sp\: UMFPACK's solve was successful.
UMFPACK V5.4.0 (May 20, 2009), Info:
matrix entry defined as: double
Int (generic integer) defined as: UF_long
BLAS library used: Fortran BLAS. size of BLAS integer: 8
MATLAB: yes.
CPU timer: POSIX times ( ) routine.
number of rows in matrix A: 4225
number of columns in matrix A: 4225
entries in matrix A: 520449
memory usage reported in: 16-byte Units
size of int: 4 bytes
size of UF_long: 8 bytes
size of pointer: 8 bytes
size of numerical entry: 8 bytes
strategy used: symmetric
ordering used: amd on A+A'
modify Q during factorization: no
prefer diagonal pivoting: yes
pivots with zero Markowitz cost: 130
submatrix S after removing zero-cost pivots:
number of "dense" rows: 0
number of "dense" columns: 0
number of empty rows: 0
number of empty columns 0
submatrix S square and diagonal preserved
pattern of square submatrix S:
number rows and columns 4095
symmetry of nonzero pattern: 1.000000
nz in S+S' (excl. diagonal): 508032
nz on diagonal of matrix S: 4095
fraction of nz on diagonal: 1.000000
AMD statistics, for strict diagonal pivoting:
est. flops for LU factorization: 1.79750e+10
est. nz in L+U (incl. diagonal): 10523205
est. largest front (# entries): 6230016
est. max nz in any column of L: 2496
number of "dense" rows/columns in S+S': 0
symbolic factorization defragmentations: 0
symbolic memory usage (Units): 1195977
symbolic memory usage (MBytes): 18.2
Symbolic size (Units): 10846
Symbolic size (MBytes): 0
symbolic factorization CPU time (sec): 0.05
symbolic factorization wallclock time(sec): 0.05
matrix scaled: yes (divided each row by sum of abs values in each row)
minimum sum (abs (rows of A)): 1.00000e+00
maximum sum (abs (rows of A)): 4.30992e+03
symbolic/numeric factorization: upper bound actual %
variable-sized part of Numeric object:
initial size (Units) 1089912 1085688 100%
peak size (Units) 15802496 11432659 72%
final size (Units) 8154726 5295040 65%
Numeric final size (Units) 8177995 5316197 65%
Numeric final size (MBytes) 124.8 81.1 65%
peak memory usage (Units) 15872229 11502392 72%
peak memory usage (MBytes) 242.2 175.5 72%
numeric factorization flops 3.65049e+10 1.79750e+10 49%
nz in L (incl diagonal) 7186015 5271844 73%
nz in U (incl diagonal) 8386756 5263908 63%
nz in L+U (incl diagonal) 15568546 10531527 68%
largest front (# entries) 12892006 6230016 48%
largest # rows in front 3149 2496 79%
largest # columns in front 4094 2496 61%
initial allocation ratio used: 0.851
# of forced updates due to frontal growth: 0
number of off-diagonal pivots: 0
nz in L (incl diagonal), if none dropped 5271844
nz in U (incl diagonal), if none dropped 5263908
number of small entries dropped 0
nonzeros on diagonal of U: 4225
min abs. value on diagonal of U: 1.00e-02
max abs. value on diagonal of U: 1.00e+00
estimate of reciprocal of condition number: 1.00e-02
indices in compressed pattern: 50236
numerical values stored in Numeric object: 10523335
numeric factorization defragmentations: 1
numeric factorization reallocations: 1
costly numeric factorization reallocations: 1
numeric factorization CPU time (sec): 9.62
numeric factorization wallclock time (sec): 2.78
numeric factorization mflops (CPU time): 1868.51
numeric factorization mflops (wallclock): 6465.84
symbolic + numeric CPU time (sec): 9.67
symbolic + numeric mflops (CPU time): 1858.84
symbolic + numeric wall clock time (sec): 2.83
symbolic + numeric mflops (wall clock): 6351.60
solve flops: 4.68820e+07
iterative refinement steps taken: 1
iterative refinement steps attempted: 1
sparse backward error omega1: 1.84e-16
sparse backward error omega2: 0.00e+00
solve CPU time (sec): 0.05
solve wall clock time (sec): 0.05
solve mflops (CPU time): 937.64
solve mflops (wall clock time): 937.64
total symbolic + numeric + solve flops: 1.80219e+10
total symbolic + numeric + solve CPU time: 9.72
total symbolic + numeric + solve mflops (CPU): 1854.11
total symbolic+numeric+solve wall clock time: 2.88
total symbolic+numeric+solve mflops(wallclock) 6257.61
justsolve =
2.8828
Solve the equations with factoring.
tic [ll,uu,pp,qq]= lu(L); factortime=toc tic u2 = qq*(uu\(ll\(pp*rhs))); uu2 = reshape(u2,N+1,N+1); factorsolve=toc
UMFPACK V5.4.0 (May 20, 2009), Control:
Matrix entry defined as: double
Int (generic integer) defined as: UF_long
0: print level: 2
1: dense row parameter: 0.2
"dense" rows have > max (16, (0.2)*16*sqrt(n_col) entries)
2: dense column parameter: 0.2
"dense" columns have > max (16, (0.2)*16*sqrt(n_row) entries)
3: pivot tolerance: 0.1
4: block size for dense matrix kernels: 32
5: strategy: 0 (auto)
6: initial allocation ratio: 0.7
7: max iterative refinement steps: 2
12: 2-by-2 pivot tolerance: 0.01
13: Q fixed during numerical factorization: 0 (auto)
14: AMD dense row/col parameter: 10
"dense" rows/columns have > max (16, (10)*sqrt(n)) entries
Only used if the AMD ordering is used.
15: diagonal pivot tolerance: 0.001
Only used if diagonal pivoting is attempted.
16: scaling: 0 (no)
17: frontal matrix allocation ratio: 0.5
18: drop tolerance: 0
19: AMD and COLAMD aggressive absorption: 1 (yes)
The following options can only be changed at compile-time:
8: BLAS library used: Fortran BLAS. size of BLAS integer: 8
9: compiled for MATLAB
10: CPU timer is POSIX times ( ) routine.
11: compiled for normal operation (debugging disabled)
computer/operating system: Linux
size of int: 4 UF_long: 8 Int: 8 pointer: 8 double: 8 Entry: 8 (in bytes)
UMFPACK V5.4.0 (May 20, 2009), Info:
matrix entry defined as: double
Int (generic integer) defined as: UF_long
BLAS library used: Fortran BLAS. size of BLAS integer: 8
MATLAB: yes.
CPU timer: POSIX times ( ) routine.
number of rows in matrix A: 4225
number of columns in matrix A: 4225
entries in matrix A: 520449
memory usage reported in: 16-byte Units
size of int: 4 bytes
size of UF_long: 8 bytes
size of pointer: 8 bytes
size of numerical entry: 8 bytes
strategy used: symmetric
ordering used: amd on A+A'
modify Q during factorization: no
prefer diagonal pivoting: yes
pivots with zero Markowitz cost: 130
submatrix S after removing zero-cost pivots:
number of "dense" rows: 0
number of "dense" columns: 0
number of empty rows: 0
number of empty columns 0
submatrix S square and diagonal preserved
pattern of square submatrix S:
number rows and columns 4095
symmetry of nonzero pattern: 1.000000
nz in S+S' (excl. diagonal): 508032
nz on diagonal of matrix S: 4095
fraction of nz on diagonal: 1.000000
AMD statistics, for strict diagonal pivoting:
est. flops for LU factorization: 1.79750e+10
est. nz in L+U (incl. diagonal): 10523205
est. largest front (# entries): 6230016
est. max nz in any column of L: 2496
number of "dense" rows/columns in S+S': 0
symbolic factorization defragmentations: 0
symbolic memory usage (Units): 1195977
symbolic memory usage (MBytes): 18.2
Symbolic size (Units): 10846
Symbolic size (MBytes): 0
symbolic factorization CPU time (sec): 0.05
symbolic factorization wallclock time(sec): 0.04
matrix scaled: no
symbolic/numeric factorization: upper bound actual %
variable-sized part of Numeric object:
initial size (Units) 1089912 1085688 100%
peak size (Units) 15802496 11432659 72%
final size (Units) 8154726 5295040 65%
Numeric final size (Units) 8177995 5314084 65%
Numeric final size (MBytes) 124.8 81.1 65%
peak memory usage (Units) 15872229 11500279 72%
peak memory usage (MBytes) 242.2 175.5 72%
numeric factorization flops 3.65049e+10 1.79750e+10 49%
nz in L (incl diagonal) 7186015 5271844 73%
nz in U (incl diagonal) 8386756 5263908 63%
nz in L+U (incl diagonal) 15568546 10531527 68%
largest front (# entries) 12892006 6230016 48%
largest # rows in front 3149 2496 79%
largest # columns in front 4094 2496 61%
initial allocation ratio used: 0.851
# of forced updates due to frontal growth: 0
number of off-diagonal pivots: 0
nz in L (incl diagonal), if none dropped 5271844
nz in U (incl diagonal), if none dropped 5263908
number of small entries dropped 0
nonzeros on diagonal of U: 4225
min abs. value on diagonal of U: 1.00e+00
max abs. value on diagonal of U: 1.87e+03
estimate of reciprocal of condition number: 5.34e-04
indices in compressed pattern: 50236
numerical values stored in Numeric object: 10523335
numeric factorization defragmentations: 1
numeric factorization reallocations: 1
costly numeric factorization reallocations: 1
numeric factorization CPU time (sec): 12.96
numeric factorization wallclock time (sec): 5.73
numeric factorization mflops (CPU time): 1386.96
numeric factorization mflops (wallclock): 3137.00
symbolic + numeric CPU time (sec): 13.01
symbolic + numeric mflops (CPU time): 1381.63
symbolic + numeric wall clock time (sec): 5.77
symbolic + numeric mflops (wall clock): 3115.26
factortime =
8.5793
sp\: bandwidth = 4195+1+0.
sp\: is A diagonal? no.
sp\: is A triangular? yes.
sp\: do a triangular solve.
sp\: bandwidth = 0+1+4067.
sp\: is A diagonal? no.
sp\: is A triangular? yes.
sp\: do a triangular solve.
factorsolve =
0.0386
The difference between the solutions.
figure(1) clf pcolor(xx,yy,uu1-uu2), shading interp, colorbar, title('Difference between solution without and with factoring')
The individual solutions.
figure(2) clf subplot(2,1,1) pcolor(xx,yy,uu1),shading interp,caxis([-1 1]) subplot(2,1,2) pcolor(xx,yy,uu2),shading interp,caxis([-1 1])
The diffusivity.
figure(3) clf pcolor(xx,yy,kappa),shading flat,colorbar, xlabel('x') ylabel('y')
Create the derivatives.
u1y=Dy*u1; u1x=Dx*u1; uy=reshape(u1y,N+1,N+1); ux=reshape(u1x,N+1,N+1);
figure(4) clf % Plot the derivatives. subplot(2,2,1) pcolor(xx,yy,ux),shading interp,caxis([-10 10]) subplot(2,2,2) pcolor(xx,yy,uy),shading interp,caxis([-10 10]) subplot(2,2,3) % Plot the derivatives at the boundaries. plot(yy(bxn1),ux(bxn1),'bo-',xx(byp1),uy(byp1),'rs-') subplot(2,2,4) plot(xx(byn1),u1(byn1),'bo-',yy(bxp1),u1(bxp1),'rs-')
Chebyshev function
CHEB compute D = differentiation matrix, x = Chebyshev grid
function [D,x] = cheb(N) if N==0, D=0; x=1; return, end x = cos(pi*(0:N)/N)'; c = [2; ones(N-1,1); 2].*(-1).^(0:N)'; X = repmat(x,1,N+1); dX = X-X'; D = (c*(1./c)')./(dX+(eye(N+1))); % off-diagonal entries D = D - diag(sum(D')); % diagonal entries end