solvefick2.m

扩散的菲克第二定律的数值求解 第一次弄

function X=SolveFick2(Matrix_C,int_IDrow,float_Cs,float_C0,r)
%int_IDrow为本次运行该函数要解决的矩阵C中行的位置（实际坐标中的列）
%float_Cs为x=0时第int_IDrow行的边界
%float_C0为x=H时第int_IDrow行的边界
%r为因子
S=size(Matrix_C);
m=S(2)-1;
n=1;
A=zeros((m+1)*(n+1));
X=zeros((m+1)*(n+1),1);
b=zeros((m+1)*(n+1),1);

%本函数中的i,j均为坐标，与矩阵中的行列的位置不同
for i=1:n
    for j=1:(m-1)
        A(Position(i,j,n,m),Position(i,j+1,n,m))=r;
        A(Position(i,j,n,m),Position(i,j,n,m))=-2*r-1;
        A(Position(i,j,n,m),Position(i,j-1,n,m))=r;
        A(Position(i,j,n,m),Position(i-1,j,n,m))=1;
    end
end

for j=0:m%t=0的情况
    A(Position(0,j,n,m),Position(0,j,n,m))=1;     
    b(Position(0,j,n,m))=Matrix_C(int_IDrow-1,j+1);
end

for i=0:n%x=0以及x=mh的情况,注意i,j均为零的情况；
    A(Position(i,0,n,m),Position(i,0,n,m))=1; 
    A(Position(i,m,n,m),Position(i,m,n,m))=1;
    b(Position(i,0,n,m))=float_Cs;
    b(Position(i,m,n,m))=float_C0;
end

X=inv(A)*b;