- 积分
- 35
- 注册时间
- 2004-10-23
- 仿真币
-
- 最后登录
- 1970-1-1
|
发表于 2006-10-15 13:13:56
|
显示全部楼层
来自 韩国
对于变扩散系数的扩散类抛物线型方程,只能采用数值方法求解。除了直接采用有限差分法和有限单元法以外,还可以采用将差分法和龙格库塔法结合起来求解。先采用差分法将二阶抛物线偏微分方程离散成为一族一阶的常微分方程族,然后才用龙格库塔法求解即可。
一个方程采用三种方法求解,就算是一个熟悉各种数值方法的过程,对于各种方法做到心中有数,不能盲目的使用软件。
对于上面的变扩散系数的抛物线型方程,下面的程序是采用差分和龙格库塔结合使用的例子:
%2006.09.19,roung kutta to solve differential equations
%xiaoyong, hy university%in this program,time unit is year ,length unit is mm,quality unit is kg
clear
clc
cs=17.7
len=200
m=50
period=50
n=100
con=zeros(1,m-1)
h=len/m
delta=period/n
han=inline('((1e-6*24*3600*365)/(1+0.98/(0.08*(1+0.29*c1)^2)))*(c0-2*c1+c2)/h^2','c0','c1','c2','h')
for i=1:n
k(1,1)=feval(han,cs,con(1),con(2),h)
for j=2:m-2
k(1,j)=feval(han,con(j-1),con(j),con(j+1),h)
end
k(1,m-1)=feval(han,con(m-2),con(m-1),0,h)
k(2,1)=feval(han,cs,con(1)+0.5*delta*k(1,1),con(2)+0.5*delta*k(1,2),h)
for j=2:m-2
k(2,j)=feval(han,con(j-1)+0.5*delta*k(1,j-1),con(j)+0.5*delta*k(1,j),con(j+1)+0.5*delta*k(1,j+1),h)
end
k(2,m-1)=feval(han,con(m-2)+0.5*delta*k(1,m-2),con(m-1)+0.5*delta*k(1,m-1),0,h)
k(3,1)=feval(han,cs,con(1)+0.5*delta*k(2,1),con(2)+0.5*delta*k(2,2),h)
for j=2:m-2
k(3,j)=feval(han,con(j-1)+0.5*delta*k(2,j-1),con(j)+0.5*delta*k(2,j),con(j+1)+0.5*delta*k(2,j+1),h)
end
k(3,m-1)=feval(han,con(m-2)+0.5*delta*k(2,m-2),con(m-1)+0.5*delta*k(2,m-1),0,h)
k(4,1)=feval(han,cs,con(1)+delta*k(3,1),con(2)+delta*k(3,2),h)
for j=2:m-2
k(4,j)=feval(han,con(j-1)+delta*k(3,j-1),con(j)+delta*k(3,j),con(j+1)+delta*k(3,j+1),h)
end
k(4,m-1)=feval(han,con(m-2)+delta*k(3,m-2),con(m-1)+delta*k(3,m-1),0,h)
for j=1:m-1
res(j)=con(j)+delta*(k(1,j)+2*k(2,j)+2*k(3,j)+k(4,j))/6
end
out(i,:)=res
%进行回带
con=res
end
axial=0:h:len
nongdu=[cs res 0]
figure
plot(axial,nongdu,'b')
gal02=[axial' nongdu']
dlmwrite('kutta100.xls',gal02,'\t')
grid
附件为有限差分法和综合法的对比,基本没有差别。
[ 本帖最后由 ilxy 于 2006-10-15 13:21 编辑 ] |
本帖子中包含更多资源
您需要 登录 才可以下载或查看,没有账号?注册
×
评分
-
1
查看全部评分
-
|