matlab有限差分法

lb107

有兄弟用过有限差分法吗？？
小弟想找几个算例看看。
谢谢！

liuzg

function cranknicolson
%This program is used to solve Partial Differential Equation by Finite Difference Formulation of Crank Nicolson Method
clc
clear
%confirm a, a has relation to time unit and quality unit,length unit
%in this program,time unit is year ,length unit is mm,quality unit is kg
hold on        %绘制的图形曲线保持，与下面运行的图形同时出现在一个图面内
a=1*1e-6*24*3600*365;    %扩散系数为常数
cs=17.7*10^9;
day=400;                        %时间，单位为年
t=200;                                %时间步长
n=day/t;                        %时间网格划分，时间步数
length=200;                %试件长度
h=20;                        %空间步长
m=length/h;                %空间网格划分
ss=a*t/(h*h);        %步长比
r=m-1;                        %空间步数
%conform coefficient matrix aaleft and bbright
aa=zeros(r,r);        %系数矩阵aa置零
bb=zeros(r,r);        %系数矩阵bb置零

%以下形成系数矩阵aa
for i=1:r
    aa(i,i)=1+ss;
end
for i=1:r-1
    aa(i,i+1)=-0.5*ss;
end
for i=2:r
    aa(i,i-1)=-0.5*ss;
end

%以下形成系数矩阵bb
for i=1:r
   bb(i,i)=1-ss;
end
for i=1:r-1
    bb(i,i+1)=0.5*ss;
end
for i=2:r
    bb(i,i-1)=0.5*ss;
end

%时间步为1时对应数值求解
u0=zeros(m-1,1);
for i=1:m-1
    u0(i)=0;
end
mul=bb*u0;
bb0=ss*0.5*(cs+cs);
bb1=ss*0.5*(0+0);
mul(1)=mul(1)+bb0;
mul(m-1)=mul(m-1)+bb1;
%trsolve is a private program used to solve three diagonal equations
u(:,1)=aa\mul;

%三对角矩阵求解
for i=1:n-1
     bb0=ss*0.5*(cs+cs);
     bb1=ss*0.5*(0+0);
    mul=bb*u(:,i);
    mul(1)=mul(1)+bb0;
    mul(m-1)=mul(m-1)+bb1;
    u(:,i+1)=aa\mul;
end

%结果可视化
u=u/10^9;
plv=[cs/10^9 u(:,n)' 0];
plh=0:h:length;
plot(plh,plv);
grid on
xlabel('Concrete Cover(mm)')
ylabel('Cf(kg/m^3)  Pore Solution')

%结果存入文件crfifty.txt
% save crfifty.txt plv plh -ascii

lb107

谢谢楼上的
楼上的能给出对应的待求的微分方程的形式吗？？

qomwfnlmm

:):(:D:P:L

qomwfnlmm

谢谢···················