为何微分方程得到的结果与文献不同

肉肉坨

方程表示为：
dP(x,t)/dt=2.3618*P(x,t)-7.81e-4*P(x,t)^3+E(x,t)+d2P(x,t)/dx2；

dE(x,t)dt=(Vref(x,t)/0.0625+35200*dVref(x,t)/dt-0.286*E(x,t)-dP(x,t)/dt)/260；

dVref(x,t)/dt=(0.1257*cos(2*pi*t)+f_int-49.2304*Vref(x,t))/1.0831e5;

其中
d2P(x,t)/dx2表示对P关于x求二次偏导

f_int表示对f=0.286*E(x,t)+dP(x,t)/dt))/260在x=(0,800)上求积分

t=(0,1)；取M=1e4；h=1/M=1e-4;
L=(0,800)；取N=200；k=800/N=4;
P(x,0)=0;
E(x,0)=0;
当x=0时，dP(x,t)/dt=0
Vref(x,0)=0;

我把f_int表示成：trapz([0:k:L],[(sigma0*E(:,1:end-1)+diff(P,1,2)/h)/epsilond])；
然后用前向差分得到的结果是：

可是文献上的结果应该是如下所示的光滑的回线？请教各位大侠那些因素可能出现这样的错误呢？

[ 本帖最后由 eight 于 2008-2-23 20:32 编辑 ]

VibrationMaster

是否步长太大了？

sigma665

两副图的坐标值不一样,他的大,不过你的怎么一圈又一圈的

肉肉坨

又重新检查了代码，发现了一些错误，可算出来的结果更加离谱了，请大家帮忙看看啊。

[ 本帖最后由 肉肉坨 于 2008-2-25 13:45 编辑 ]

肉肉坨

1. clear;
   
2. Cref=22; %reference capacitor
   
3. Aferro=6.25e-4; %capacitor areas
   
4. Rref=100; %reference resistence
   
5. G=1; %gradiet constant
   
6. V0max=0.02; %V0 max value
   
7. epsilond=260; %d layer
   
8. epsilonf=260; %f layer
   
9. sigma0=2.86e-7; %conductivity
   
10. alphaf=-2.3618; %2 order constant in flayer
    
11. beltaf=7.81e-4; %4 order constant in f layer
    
12. alphad=-236.18; %2 order constant in d layer
    
13. beltad=7.81e-2; %4 order constant in d layer
    
14. nu=0.1; %width ratio
    
15. L=800; %width
    
16. Ld=L*nu; %d layer width 8e-6
    
17. Lf=L*(1-nu); %f layer width
    
18. N=800;
    
19. k=L/N; %△x:
    
20. epsilon=ones(N+1,1);
    
21. epsilon(1:nu*N)=epsilond;
    
22. epsilon(nu*N+1:end)=epsilonf;
    
23. T=1;
    
24. M=1e3;
    
25. h=T/M; %△t:
    
26. 
    
27. %-------------initial boundary condition----------%
    
28. Vref=zeros(N+1,M+1);
    
29. fVref_integrate=zeros(N+1,M+1);
    
30. E=zeros(N+1,M+1);
    
31. dV0=ones(N+1,M+1);
    
32. P=zeros(N+1,M+1);
    
33. fP=zeros(N+1,M+1);
    
34. P(:,1)=0;
    
35. P(:,2)=P(:,1);
    
36. Eave=zeros(1,M+1);
    
37. Pave=zeros(1,M+1);
    
38. Vave=zeros(1,M+1);
    
39. 
    
40. %------------forward eular method II-----------------%
    
41. for j=1:M %time difference
    
42. dV0(:,j)=repmat(cos(2*pi*j*h/1000),N+1,1);
    
43. %fVref_integrate(:,j)=repmat(trapz([0:k:L],[(sigma0*E(:,j)+diff(P(:,j:j+1),1,2)/h)./epsilon]),N+1,1);
    
44. fVref_integrate(:,2:end)=repmat(trapz([0:k:L],[(sigma0*E(:,1:end-1)+diff(P,1,2)/h)/epsilond]),N+1,1); %求积分
    
45. Vref(:,j+1)=h*(V0max*2*pi*dV0(:,j)+fVref_integrate(:,j)-(Ld/epsilond+Lf/epsilonf)*Vref(:,j)/Rref/Aferro)/(1+Cref/Aferro*(Ld/epsilond+Lf/epsilonf))+Vref(:,j);
    
46. E(:,j+1)=h*(Vref(:,j)/Rref/Aferro+Cref/Aferro*(Vref(:,j+1)-Vref(:,j))/h-sigma0*E(:,j)-(P(:,j+1)-P(:,j))/h)/epsilond+E(:,j);
    
47. for i=2:N %spatial difference
    
48. fP(i,j)=-alphaf*P(i,j)-beltaf*P(i,j)^3+E(i,j)+G*(P(i+1,j)+P(i-1,j)-2*P(i,j))/(k^2); %中心差分格式
    
49. end
    
50. P(:,j+1)=h*fP(:,j)+P(:,j);
    
51. end
    
52. Eave=mean(E);
    
53. Pave=mean(P);
    
54. Vave=mean(Vref);
    
55. plot(Eave);
    
56. ylabel('Eave');
    
57. figure;
    
58. plot(Pave);
    
59. ylabel('Pave')
    
60. figure;
    
61. plot(Vave);
    
62. ylabel('Vave')
    
63. figure;
    
64. plot(Eave,Pave);
    
65. xlabel('Eave');ylabel('Pave');
    
66. figure;
    
67. plot(Eave,Vave);
    
68. xlabel('Eave');ylabel('Vave');

复制代码

我试图用ode45等MATLAB中的语句来写一阶的常微分的方程组，但是因为方程的右边又有微分，又有积分，都不知道怎么写odefun函数了，所以就自己写了前向差分的欧拉法程序。

[ 本帖最后由 肉肉坨 于 2008-2-25 13:51 编辑 ]

无水1324

直接运行之后，图跟你的图右很大的差异

肉肉坨

你的图是不是这样的呢？

我再把式子列出来，是这样的：

其中的未知数是P，E，Vref。其余的所有参数均已知。epsilon是分段函数，在x=(0,Ld)时，epsilon=epsilond。在x=(Ld,L)时，epsilon=epsilonf。Jc(x,t)=sigma0*E(x,t).v0=V0max*sin(ωt).
文献中给出的差分参数分别是N=200，M=1e6；我的计算机算不了那么多时间步。我把M减小为1000，并且适当增加了v0的频率。
得到的图形还不是文献上的图的样子。我想问的是，这个算法有问题吗?到底步长要怎么取才好呢

[ 本帖最后由 肉肉坨 于 2008-2-27 16:06 编辑 ]