1
$\begingroup$

for a convection-diffusion equation with Dirichlet boundary conditions as follows: $$-u''+qu'=f$$ Using centered difference for $u''$ and $u'$, we get a linear system $$Ax=b$$where matrix $A$ is nonsymetric. Let $H = (A+A')/2,S = (A-A')/2$, we know that matrix $H$ is symmetric positive definite. Given a constant $\alpha>0$, we can get two linear large sparse systems, and $I$ is a identity matrix. $$(\alpha I+H)x_{k+1/2} = b_1$$and $$(\alpha I+S)x_{k+1} = b_2$$

My question is in 3D case, the above two systems are so large and sparse, how to use sin and modified sin transfomation to solve the 2 systems quickly, or some other fast methods, because in matlab, the backslash command '\' tend to be slow for this very large sparse in 3D. Can someone give the matlab codes, thanks very much.

Below is my matlab codes: for n=64 equidistance so slow.

clc;clear;close all
n=64;q=1;
h=1/(n+1);
I = speye(n);
r = q*h/2;
t1 = 6; t2 = -1-r; t3 = -1+r;
Tx = spdiags(ones(n,1).*[t2 t1 t3],[-1 0 1],n,n);
Ty = spdiags(ones(n,1).*[t2 0 t3],[-1 0 1],n,n);
Tz = spdiags(ones(n,1).*[t2 0 t3],[-1 0 1],n,n);
A = kron(Tx,kron(I,I))+kron(I,kron(Ty,I))+kron(I,kron(I,Tz));
H = (A+A')/2;   S = (A-A')/2;
lambda_maxH = eigs(H,1,'largestabs');
lambda_minH = eigs(H,1,'smallestabs');
alpha = sqrt(lambda_maxH*lambda_minH);
N = length(A);
I = speye(N);
b = A*ones(N,1);

%   this 2 systems are solved so slow
x1 = (alpha*I+H)\b;
x2 = (alpha*I+S)\b;

$\endgroup$
3
  • 1
    $\begingroup$ I would look at solving these equations iteratively (say with MINRES/GMRES). You will probably need a good preconditioner. There is a ton of literature on preconditioners for these multi-level Toeplitz systems, but I don't think there a silver bullet. A good introduction is: pdfs.semanticscholar.org/5138/… $\endgroup$ Oct 28 '19 at 8:04
  • $\begingroup$ thanks Prof for your reply, I know that we would better choose iterative methods. But if we want to solve these two large shifted systems w.r.t. (\alpha I +H) and (\alpha I +S) accurately, can the fast fourier method work? and how it do that? thanks very much. $\endgroup$
    – sunshine
    Oct 28 '19 at 8:23
  • 1
    $\begingroup$ I'm not quite sure what you mean by the Fourier method. The basic preconditioner for a Toeplitz system is the Strang preconditioner, in which the Toeplitz matrix is approximated by a circulant matrix. The FFT is used for solving the circulant system, as it is diagonal in the Fourier basis. It is also used for computing the matrix-vector product fast when the Toeplitz matrix is dense (in your case this isn't a problem). The sine transformation is a related idea; it is explained in the paper I linked to. $\endgroup$ Oct 28 '19 at 9:34

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.