[求助]偏微分方程边界条件随时间变化的问题

karie

我有一个偏微分方程边界条件是随时间变化的,请问怎么才能解,我试了好多方法!!谢谢指导.

1. clear
   
2. kongjian=1;
   
3. shijian=50000;
   
4. t=1000;
   
5. h=0.02;
   
6. n=(kongjian/h)+1;
   
7. m=(shijian/t)+1;
   
8. cs=2;
   
9. k=0.001;
   
10. %k=1;
    
11. k1=0.000001
    
12. u1=zeros(n,m);
    
13. %a xiao c de kuo san xi shu
    
14. a=0.0000001;
    
15. r=1/(t*k/2+1);
    
16. s=a*t/h^2;
    
17. ug=20;
    
18. lamud=5e-008;              
    
19. %bian jie tiao jian
    
20. %u1(1,1:m)=20;
    
21. % chu shi tiao jian
    
22. %con=0
    
23. u1(1:n,1:m)=1;
    
24. con=4
    
25. 
    
26. for j=2:m
    
27.       u1(1,j)=r*(a*t/h^2*(u1(2,j-1)-2*u1(1,j-1)+u1(1,j-1)+h*lamud/a*(ug-u1(1,j-1)))+(1-t*k/2)*u1(1,j-1)+t*k*cs)
    
28. 
    
29.     for i=2:n-1
    
30.        % if h*(i)<=(k1*1000*(j-1))^0.5
    
31.             if i<con-1                              
    
32.            u1(i,j)=r*(a*t/(h^2)*(u1(i+1,j-1)-2*u1(i,j-1)+u1(i-1,j-1))+(1-t*k/2)*u1(i,j-1)+t*k*cs)
    
33.              else
    
34.                  if con-1==i
    
35.                         u1(i,j)=r*(a*t/h^2*(u1(i,j-1)-2*u1(i,j-1)+u1(i,j-1)+h*lamud/a*(u1(i-1,j-1)-u1(i,j-1)))+(1-t*k/2)*u1(i,j-1)+t*k*cs)
    
36.           % continue
    
37.                 else  
    
38.                         u1(i,j)=r*(a*t/(h^2)*(u1(i+1,j-1)-2*u1(i,j-1)+u1(i-1,j-1))+(1-t*k/2)*u1(i,j-1)+t*k*cs)   
    
39.            %continue
    
40.      end
    
41. end
    
42.      if h*(i)>(k1*t*(j-1))^0.5
    
43.             con=i
    
44.             break
    
45.         else
    
46.             %break
    
47.           continue
    
48.       end
    
49.     break
    
50.     u1(n,j)=u1(n-1,j)
    
51. end
    
52. u1(n,j)=u1(n-1,j)
    
53. end
    
54. w=linspace(0,shijian,m);
    
55. w2=u1(1,:)
    
56. w3=u1(2,:)
    
57. w4=u1(4,:)
    
58. w5=u1(6,:)
    
59. w6=u1(8,:)
    
60. w8=u1(10,:)
    
61. %plot(w,w2,'b')
    
62. ww=linspace(0,kongjian,n);
    
63. ww2=u1(:,40)
    
64. ww3=u1(:,10)
    
65. ww4=u1(:,2)
    
66. 
    
67. %plot(ww,ww2,'r')
    
68. subplot(2,1,1),plot(w,w2,'b',w,w3,'r',w,w4,'y',w,w5,'c',w,w6,'k',w,w8,'g')
    
69. subplot(2,1,2),plot(ww,ww2,'r',ww,ww3,'b',ww,ww4,'y')
    
70. 
    
71. %fevel(pde,n-haha,m-hehe)

复制代码

happy

什么地方有问题？不是可以运行吗？

karie

是可以运行,但是结果不和逻辑,因为是个反应扩散方程.所以算出来的答案u1横向应该增大,纵向减小.但是因为边界条件的关系总是有的数据不对.不知道边界条件随时间变化,编程的时候是不是有这方面的好方法.谢谢教授!这里好像不能加附件,那我在发表一个新帖子,方程和边界条件放在那里.

happy

你这个问题好像没这么简单
看了一下你的出程序，你用的是差分吧，不知道你的计算网格是怎么给出来的呢？

karie

我用的显示差分格式,0<x<1cm,时间t=50000秒,空间步长我用的是0.02,时间步长1000,空间步长和时间步长的取值是我参照别人渗碳方程取的,因为时间太长,所以时间步长取值不可能很小.谢谢教授指导!

happy

能否先介绍一下空间坐标你是如何划分的

karie

<P>因为是个一维问题,所以只有x一个变量,0&lt;x&lt;1,分成50份,每份大小0.02,程序中我用x来表示.在t=0时,边界设在x=0处,然后边界随时间变化遵照我word文档中的关系逐渐向前走.谢谢教授.</P>

karie

<P>请帮忙,万分感谢!!,急!!</P>

happy

这个能否仔细说明一下
换句话说就是你的差分方程式如何构造的

karie

<P>谢谢happy来看,前两天网络有问题,没上来.我的差分是用的就是普通的显示格式,我把算法写贴楼主层,放在附件里,谢谢!</P>

happy

程序上应该没有他大问题，除非你把符号之类的输错了那就另当别论了

总觉得如果出问题应该在算法上，或者离散上，不知道你算没算过特征线？

karie

<P>非常感谢,我没有用过特征线,不会用.这个方程的边界条件不固定,给我带来了很大的困难,我总是觉得这个程序算出的数据有问题.这种算变化边界条件的差分是不是还有别的方法,有没有固定的套路可寻那,或是边界条件怎么加比较好那,请教happy,多谢指导!!</P>

happy

这方面我做的不多，上次也碰到类似变边界条件的问题，后来请教了一个研究流体的老师，他说我问的那个问题他正要立项做作，后来就放弃了

问题的关键还是差分的构造问题，分析特征线是必须的，只有这样才有可能构造出合理的差分格式，否则结果可能就有可能相差很大，或者可以说面目全非吧

由于这方面几乎没什么经验，不能给你提供什么帮助，请见谅

karie

<P>谢谢happy这么长时间的帮助.非常感激.:)</P>

songandyou

   可以用有限差分法解吧