matlab有限差分法
有兄弟用过有限差分法吗？？小弟想找几个算例看看。
谢谢！
参考参考吧
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=;
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 谢谢楼上的
楼上的能给出对应的待求的微分方程的形式吗？？ :):(:D:P:L 谢谢···················
页:
[1]