Help! Why is my code so slow/ does this even work properly? Using a crank-Nicholson scheme to solve a PDE.

2 次查看（过去 30 天）
Cathleen
Cathleen 2012-4-9
Hello,
I am using a numerical scheme (crank-Nicholson scheme) to solve a PDE as part of an assignment.
This is an diffusion-advection-reaction problem.
The equation in question is the following.
dC/dt= D*[d^2*C/dx^2]- ad*[dC/dx]- R(x)
The first term and second terms have constant coefficients, while the reaction term is a function of x. Part of using the crank Nicholson scheme is inverting a matrix.As you can see, with a great number of points this can be a problem (when using the \ command). I tried using the conjugate gradient code provided by mathworks, but either 1) I am not using it correctly or 2) the way my code works makes everything slow. The code is below. Any tips or suggestions will be appreciated!
~MsCaptainFlan **************************************************************
function [pp ] = CrankNick3
%Solves the advection-diffusion-reaction PDE using the crank-nicholson
%scheme
% instead of using the command AA\BB, the conjugate graduent method is used
% to speed up computation
%the conjugate gradient code was obtained from MathWorks
xl=0;xr=1;%boundaries in space
yb=0;yt=1; %boundaries in time
M=50; %number of points in space
N=100; %number of points in time
dx=(xr-xl)/M;dt=(yt-yb)/N; % step sizes
m=M-1; n=N;
x=xl:dx:xr;tt=yb:dt:yt;
xx=x(1:m);
G0=10; %uM
%coefficients
global D W k b tol
D=8e-8; %diffusion
W=-1*.1e-2;%advection
tol=1e-5;
%Rxn should be a m*1 vector. The rxn coefficient of the PDE varies in space
k=.15; %Rcox from Burdige
b=.1; %the beta reaction term
Rxn=k*exp(-b*xx');
%Coefficients used to generate matrix AA and BB
global alpha beta gama
alpha=dt*D/(2*dx);
beta=-dt*W/(4*dx);
gama=-dt*Rxn./2;
global c1 c2 c3 c4 c5 c6
c1=-alpha-beta;
c2=1+2*alpha-gama(1:m);
c3=-alpha+beta;
c4=alpha+beta;
c5=1-2*alpha+gama(1:m);
% c5=gama(1:m);
c6=alpha-beta;
%Generate AA matrix
% define tridiagonal matrix AA
AA=diag(c2.*ones(m,1))+diag(c3.*ones(m-1,1),1);
AA=AA+diag(c1*ones(m-1,1),-1);
%the next few lines define the boundary conditions
AA(1,:)=0; %the first row of any colum must be zero
AA(1,1)=1; %the first row, first column value must be one
AA(m,:)=0; %the last row, any colum must be zero
AA(m,m)=1; %the last row, last column must be one
%form initial matrix B0 with some concentration profile
% B0=G0;
Bb=zeros(m,1); %this forms the initial b matix
Bb(1,1)=(G0); %this defines the initial conditions
w0=AA\Bb;
BB=zeros(m,1);
BB=Bb; %for the first iteration of BB, let it be initial condition vector Bb
w=w0; %for the first iteration of w, let it be initial condition vector w0
tic
for j=1:n,
for p=2:m,
BB(p)=(c4*w(p))+(c5(p).*w(p))-(c6*w(p));
BB;
gcf;
xp=x';
fignodiff=plot(xp(1:49,1),w);hold on;
xlabel('depth');ylabel('c');
drawnow;
end
w=conjgrad(AA,BB,tol);
end
toc
end
function w = conjgrad(AA,BB,tol)
% CONJGRAD Conjugate Gradient Method.
% X = CONJGRAD(A,B) attemps to solve the system of linear equations A*X=B
% for X. The N-by-N coefficient matrix A must be symmetric and the right
% hand side column vector B must have length N.
%
% X = CONJGRAD(A,B,TOL) specifies the tolerance of the method. The
% default is 1e-10.
%
% Example:
%{
n = 6000;
m = 8000;
A = randn(n,m);
A = A * A';
b = randn(n,1);
tic, x = conjgrad(A,b); toc
norm(A*x-b)
%}
% By Yi Cao at Cranfield University, 18 December 2008.
%
if nargin<3
tol=1e-10;
end
r = BB + AA*BB;
y = -r;
z = AA*y;
s = y'*z;
t = (r'*y)/s;
w = -BB + t*y;
for k = 1:numel(BB);
r = r - t*z;
if( norm(r) < tol )
return;
end
B = (r'*z)/s;
y = -r + B*y;
z = AA*y;
s = y'*z;
t = (r'*y)/s;
w = w + t*y;
end
end

请先登录，再回答此问题。

回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Partial Differential Equation Toolbox 的更多信息

标签

产品

Community Treasure Hunt

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

Start Hunting!

Translated by