常偏微分方程解高手请参与讨论!

bainhome

终于连练带说地，算是有点儿PDE了^_^
PDETOOL是个很好的工具，作为一种教学和自学的手段，对于抽象概念的理解比较有好处。PDE求解由于没有自己编过程序，不便发言。而一般玩儿偏微分方程的都喜欢用femlab，我装了一回，自己的机器太慢，又卸掉了。
anyway，欢迎继续讨论！

[ 本帖最后由 bainhome 于 2006-6-3 16:07 编辑 ]

ilxy

学习matlab的过程基本上是看着simwe上的帖子学过来的，现在针对一个抛物线偏微分方程简单的例子，说明一下解析法，有限差分法以及有限单元法的基本过程。也希望每一个网友都能象bainhome一样热心。
part1：
1.参考文献< a study of the effect of chloride binding on service life predictions>(<cement and concrete research>30(2000),1215-1223, B.Martin-Perez,etc)
2。差分法采用.CRANK-NICOLSON scheme，有限单元法的过程时间和空间都采用伽辽金方法离散。
3.初始条件：
  时间t＝0时，cf＝0
  边界条件：
    x＝0时，cs＝17.7
   x＝L时。cs＝0
  计算试件长度L=200
4.方程的形式： dc/dt=D(d2c/dx2),其中D为扩散系数
5。该方程存在解析解，解析解的求解过程详见各种参考书。
6。差分法源程序：
function crank
%this program is made by xiaoyong and can be used to solve partia
%differential equation by finite difference formulation 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
a=1*10^(-6)*24*3600*365
cs=17.7*10^9
day=50
t=2
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)
bb=zeros(r,r)
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

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
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
xlabel('concrete cover(mm)')
ylabel('cf(kg/m3 pore solution')
save crfifty.txt plv plh -ascii
7.有限单元发源程序
%galerkin method,diffusion filed,think carefully
%hanyang university,erc lab,2006.06.06,wang xiaoyong
%in this program,time unit is year ,length unit is mm,quality unit is kg
clear
clc
a=1*10^(-6)*24*3600*365
cs=17.7
length=200
elnum=10
ellen=length/elnum
nonum=elnum+1
nofreedom=1
freedom=nonum*nofreedom
he=[1/ellen,-1/ellen;
    -1/ellen,1/ellen]
ce=[ellen/3,ellen/6;
    ellen/6,ellen/3]*1/a
%矩阵集成
h=zeros(freedom,freedom)
c=zeros(freedom,freedom)
%生成单元定位向量
dir(1,:)=1:2
for i=2:elnum
    dir(i,:)=dir(i-1,:)+1
end
%按照单元定位向量进行定位,实行分块定位的方法
h(dir(1,:),dir(1,:))=he
for i=2:elnum
    h(dir(i,:),dir(i,:))=h(dir(i,:),dir(i,:))+he
end
c(dir(1,:),dir(1,:))=ce
for i=2:elnum
    c(dir(i,:),dir(i,:))=c(dir(i,:),dir(i,:))+ce
end
%对于边界条件进行处理，对于h和c进行降维,
ff=zeros(freedom-2,2)
for i=1:freedom-2
ff(i,1)=h(i+1,1)
ff(i,2)=h(i+1,freedom)
end
%构造右侧矩阵ff
tt=[cs;0]
ffright=-1*ff*tt

h=h(2:freedom-1,2:freedom-1)
c=c(2:freedom-1,2:freedom-1)
alafa=2/3
day=50
deltime=2
%the initial condation is: cs 000
x0=zeros(1,elnum-1)'
aaleft=c+h*alafa*deltime
bbright=c-h*(1-alafa)*deltime

for i=1:day/deltime
    x1=inv(aaleft)*(bbright*x0+ffright*deltime)
    x0=x1
end
%最后结果绘图
x2=[cs;x1;0]
x2=x2'
axial=0:ellen:length
plot(axial,x2,'b')
grid
8。差分法和有限元法计算结果的对比，比较结果为50年时的氯离子浓度：
有限单元发结果：x2 =

   17.7000   12.8695    8.5786    5.1961    2.8350    1.3806    0.5938    0.2227    0.0715    0.0185         0
差分法结果：plv =

   17.7000   12.7532    8.4118    5.0598    2.7727    1.3859    0.6334    0.2655    0.1016    0.0330         0
9。这是一个简单的例子，有限元法和差分法太大的区别。当偏微分方程组耦合性很复杂时，差分法很难进行，推荐使用有限单元法。
10。附图：差分法和解析法计算结果对比：

[ 本帖最后由 ilxy 于 2006-6-8 10:43 编辑 ]

bainhome

我发现ilxy钻研精神非常值得称道，不愧是okok的版主，很多问题往往独立参照文献自己实现，虽然现在的MATLAB仍然有很多还可以提高的，但是直接联想专业问题进行实现的思考方法我感觉已经望尘莫及^_^，因为我们是同行，所以惭愧一下hoho...
另外提个小意见，其中有许多循环直觉上认为可以矢量优化一下，比如：

1. for i=1:r
   
2.     aa(i,i)=1+ss
   
3. end
   
4. for i=1:r-1
   
5.     aa(i,i+1)=-0.5*ss
   
6. end
   
7. for i=2:r
   
8.     aa(i,i-1)=-0.5*ss
   
9. end
   
10. 
    
11. for i=1:r
    
12.    bb(i,i)=1-ss
    
13. end
    
14. for i=1:r-1
    
15.     bb(i,i+1)=0.5*ss
    
16. end
    
17. for i=2:r
    
18.     bb(i,i-1)=0.5*ss
    
19. end
    
20. u0=zeros(m-1,1)
    
21. for i=1:m-1
    
22.     u0(i)=0
    
23. end
    

复制代码

ilxy

想了一下，采用分块矩阵，可以这么改了：
bb=diag(ones(r,1)*ss+1)
aa=diag(ones(r-1,1)*ss*0.5)
aa1=[zeros(1,r);aa,zeros(r-1,1)]
cc1=[zeros(r-1,1),aa;zeros(1,r)]
aaleft=bb+aa1+cc1
这样就把三对角矩阵构造出来了，呵呵

WaitingForMe

求解一些复杂的方程或者偏微分方程往往不是短时间内可以解决的。这个是非常实际的问题。一些参考书上的例题，实在太简单，照着学习了一遍估计也不会有太大的实质性进步。很多来源于实际问题的方程求解，最优化问题，常微分偏微分问题，具有非常复杂的工程背景，如果没有花一定时间去做基础研究，得到的解很有可能是不合理的解。而大多数人来这里，来也匆匆，去也匆匆，我想不太可能一眼望去问题就立马解决。我想这个是导致这类问题难以得到解决的真正原因。

      我想要想真正的觉得这个问题，单单凭大家的一腔热情，即使像bainhome这样热心的大侠，我想也难以支持太久。毕竟，一个人能力越强，往往自己肩上的工作也越重，难以将自己的精力集中到这样一些与自己毫不相干而且难度较高的问题上面来。我想，能够做的，也就是碰到自己恰好在做的问题，给出一些自己的经验或者教训，这样的人已经很不错了。

       解决问题的一个途径是，寻找与自己的课题类似的人互相探讨交流，一起解决问题。毕竟现在的学科分的太细，每个学科又各自有自己的一段历史的积淀，要想让学别的学科的人来帮助解决这样那难度上更是增大了很多。一些基础问题很多人都知道，也很愿意回答，但是对一些难度高专业性强的问题，不要期盼太多，还是自己静下心来好好研究。

      另一个解决的办法，可以有偿的向论坛里面的高人请求帮助。原因是，你没有理由不劳而获。或者说，以不劳而获为耻。面对一些高难度的问题，我看很多人往往一边伸手要别人帮助，一边抱怨世态炎凉。我想这样的人跟乞丐没有什么区别。每个人都要生活，你有理由让别人放下自己的饭碗来解决这样一些问题吗？有句话叫重金之下必有勇夫。

        以上是我个人的一些愚见。

bainhome

原帖由 ilxy 于 2006-6-14 10:26 发表
想了一下，采用分块矩阵，可以这么改了：
bb=diag(ones(r,1)*ss+1)
aa=diag(ones(r-1,1)*ss*0.5)
aa1=
cc1=
aaleft=bb+aa1+cc1
这样就把三对角矩阵构造出来了，呵呵

========================================================================
构造三对角阵再提供一种方法，其实就是换换diag的次序^_^

1. a=[1:6];b=[6:-1:1];c=[2:8];
   
2. d=diag(a,1)+diag(b,-1)+diag(c)

复制代码

不知道是不是ilxy想实现的三对角阵构造更理想方式:)
========================================================================
to waitingforme：一次写这么多，而且没多少让人看着不顺眼的话，太不容易了，对你...:lol

[ 本帖最后由 bainhome 于 2006-6-15 01:48 编辑 ]

ilxy

对于变扩散系数的扩散类抛物线型方程，只能采用数值方法求解。除了直接采用有限差分法和有限单元法以外，还可以采用将差分法和龙格库塔法结合起来求解。先采用差分法将二阶抛物线偏微分方程离散成为一族一阶的常微分方程族，然后才用龙格库塔法求解即可。
一个方程采用三种方法求解，就算是一个熟悉各种数值方法的过程，对于各种方法做到心中有数，不能盲目的使用软件。
对于上面的变扩散系数的抛物线型方程，下面的程序是采用差分和龙格库塔结合使用的例子：
%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 编辑 ]

ilxy

一个ansys计算扩散问题的实例
http://www.simwe.com/forum/thread-806873-1-1.html
有限元解法在ansys中的具体实现

wandererdong

I=E-A.dI/dt    I为要求解的值（N个变量），dI/dt为I的微分。这里E和A为构造出的 N X 1和N X N矩阵。有人采用离散后求解，并将delta t作为计算步长，但是没有搞明白其中的求解方法，而且文献有些前后矛盾。哪位高手给解一下？

face05

这么好的帖子应该让置顶啊，搜了一晚上，终于把他找出来了，受益颇多

axi.001

不过实际问题一般都是大量的方程组的求解啊，而且有大量的参数偶合，相当麻烦，
该讨论一下方程组的问题

coblan

我来发发我的夭折了的论文上的关于方程组的微分方程解法。我的论文是利用有限元的基础知识就汽车的二分之一模型离散化。利用拉格朗日得出的动力学方程为：
MX''+CX'+KX=R (其中，M是质量矩阵，C是阻尼矩阵，K是刚度矩阵，R是输入向量)
M,C,K可以根据耗散能量和势能动能，计算得出。
由于微分方程存在刚性问题，所以选用了隐式R-K法

coblan

以上是将R-K法首先化简为一个等式，这样在编写程序时，可以大大的简化代码，也可以充分的利用matlab矩阵运算能力。由于汽车模型的参数非常多，这里用了很多全局变量。
function[tt,y,v,a]=rk_4_yingshi_ding(n,h)
global kf0 kr0 kf kr z0f z0r z1f z1r
global M1 dem check check2
global cf cr

I1=eye(dem);
y=zeros(dem,n);tt=zeros(1,n);v=zeros(dem,n);a=zeros(dem,n);
t=0;
I=eye(2*dem);
X=zeros(2*dem,1);        
for i=1:n
C=CC(cf,cr);
kf0=604685.4;
kr0=985965;
% if X(dem+z1f)-z0f>6.6473e3/604685.4 %前轮离地判断
%     kf0=6.6473e3/(X(dem+z1f)-z0f);
%   %  check(i)=(X(dem+z1f)-z0f)-6.6473e3/604685.4;
% end
% if X(dem+z1r)-z0r>5.1735e3/985965  %后轮
%     kr0=5.1735e3/(X(dem+z1r)-z0r);
%    % check2(i)=(X(dem+z1r)-z0r)-5.1735e3/985965;
% end
%  check(i)=6.6473e3-kf0*(X(dem+z1f)-z0f);
% check2(i)=5.1735e3-kr0*(X(dem+z1r)-z0r);

K=KK(kf0,kf,kr0,kr);
Z=[-M1*C,-M1*K;I1,0*I1];
ZZ=[I-h/4*Z,(sqrt(3)/6-1/4)*h*Z;-(1/4+sqrt(3)/6)*h*Z,I-h/4*Z]^-1;

t=t+h;
R=RR(t);
Rm=[M1*R;zeros(dem,1)];

    X=X+h/2*[I,I]*ZZ*([Z*X;Z*X]+[Rm;Rm]);
%==========下面是提取该循环的输出量==========
y(:,i)=X(dem+1:2*dem);
v(:,i)=X(1:dem);
tt(1,i)=t;
a(:,i)=[-M1*C,-M1*K]*X+M1*R;
end

end

coblan

我代码的意思是求解微分方程组式，先化解欧拉，R-K，adams，gear,等等公式，再利用循环计算微分方程。如果非线性的问题，可以在循环内形成M C K矩阵，如果是线性，或者M C K的情况比较少，可以在循环外形成，这样遇到很多自由度的模型时，可以极大的提高运算速度。

bainhome

抱歉，很长时间才回复，前面一直漏看了这个帖子，我觉得您给出的推导，基本思想是状态空间的思想，在解决常系数微分方程的问题时，这的确是行之有效的。
赞一个！

liuchanghong

看完整个帖子，觉着各位都很强大。弱弱地发表一下体会。我是做混沌时使用MATLAB编程的，但是当时感觉效果不满意，后来请教前辈师兄，教给我使用他自编的龙格库塔的fortran程序，做完课题后就再也没有用MATLAB解偏微分方程。

I小浣熊

学习了~~~高手很多啊~~

lovezhiyu

大家好 我是初学者 看到楼上各位给出的解很佩服 不知道有没有对符号常数偏微分方程有研究的 希望不吝赐教 我遇到了个难题

伤痕累累

WaitingForMe 发表于 2006-6-14 13:35
求解一些复杂的方程或者偏微分方程往往不是短时间内可以解决的。这个是非常实际的问题。一些参考书上的例题 ...

i do agree with you. i'm trouble in it  for a long time for nonliner PDE.and i really want someelites to conquer it altogether !

伤痕累累

WaitingForMe 发表于 2006-6-14 13:35
求解一些复杂的方程或者偏微分方程往往不是短时间内可以解决的。这个是非常实际的问题。一些参考书上的例题 ...

i do agree with you. i'm trouble in it  for a long time for nonliner PDE.and i really want someelites to conquer it altogether !