Help! Why is my code so slow/ does this even work properly? Using a crank-Nicholson scheme to solve a PDE.
2 次查看(过去 30 天)
显示 更早的评论
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 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Boundary Conditions 的更多信息
产品
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!