- 积分
- 0
- 注册时间
- 2006-6-14
- 仿真币
-
- 最后登录
- 1970-1-1
|
发表于 2006-8-15 13:56:19
|
显示全部楼层
来自 江苏无锡
参考参考吧
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 |
|