为何微分方程得到的结果与文献不同
方程表示为：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(,[(sigma0*E(:,1:end-1)+diff(P,1,2)/h)/epsilond])；
然后用前向差分得到的结果是：
可是文献上的结果应该是如下所示的光滑的回线？请教各位大侠那些因素可能出现这样的错误呢？
[ 本帖最后由 eight 于 2008-2-23 20:32 编辑 ] 是否步长太大了？
回复 楼主 的帖子
两副图的坐标值不一样,他的大,不过你的怎么一圈又一圈的 又重新检查了代码，发现了一些错误，可算出来的结果更加离谱了，请大家帮忙看看啊。[ 本帖最后由 肉肉坨 于 2008-2-25 13:45 编辑 ] clear;
Cref=22; %reference capacitor
Aferro=6.25e-4; %capacitor areas
Rref=100; %reference resistence
G=1; %gradiet constant
V0max=0.02; %V0 max value
epsilond=260; %d layer
epsilonf=260; %f layer
sigma0=2.86e-7; %conductivity
alphaf=-2.3618; %2 order constant in flayer
beltaf=7.81e-4; %4 order constant in f layer
alphad=-236.18; %2 order constant in d layer
beltad=7.81e-2; %4 order constant in d layer
nu=0.1; %width ratio
L=800; %width
Ld=L*nu; %d layer width 8e-6
Lf=L*(1-nu); %f layer width
N=800;
k=L/N; %△x:
epsilon=ones(N+1,1);
epsilon(1:nu*N)=epsilond;
epsilon(nu*N+1:end)=epsilonf;
T=1;
M=1e3;
h=T/M; %△t:
%-------------initial boundary condition----------%
Vref=zeros(N+1,M+1);
fVref_integrate=zeros(N+1,M+1);
E=zeros(N+1,M+1);
dV0=ones(N+1,M+1);
P=zeros(N+1,M+1);
fP=zeros(N+1,M+1);
P(:,1)=0;
P(:,2)=P(:,1);
Eave=zeros(1,M+1);
Pave=zeros(1,M+1);
Vave=zeros(1,M+1);
%------------forward eular method II-----------------%
for j=1:M %time difference
dV0(:,j)=repmat(cos(2*pi*j*h/1000),N+1,1);
%fVref_integrate(:,j)=repmat(trapz(,[(sigma0*E(:,j)+diff(P(:,j:j+1),1,2)/h)./epsilon]),N+1,1);
fVref_integrate(:,2:end)=repmat(trapz(,[(sigma0*E(:,1:end-1)+diff(P,1,2)/h)/epsilond]),N+1,1); %求积分
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);
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);
for i=2:N %spatial difference
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); %中心差分格式
end
P(:,j+1)=h*fP(:,j)+P(:,j);
end
Eave=mean(E);
Pave=mean(P);
Vave=mean(Vref);
plot(Eave);
ylabel('Eave');
figure;
plot(Pave);
ylabel('Pave')
figure;
plot(Vave);
ylabel('Vave')
figure;
plot(Eave,Pave);
xlabel('Eave');ylabel('Pave');
figure;
plot(Eave,Vave);
xlabel('Eave');ylabel('Vave');
我试图用ode45等MATLAB中的语句来写一阶的常微分的方程组，但是因为方程的右边又有微分，又有积分，都不知道怎么写odefun函数了，所以就自己写了前向差分的欧拉法程序。
[ 本帖最后由 肉肉坨 于 2008-2-25 13:51 编辑 ]
回复 5楼 的帖子
直接运行之后，图跟你的图右很大的差异 你的图是不是这样的呢？我再把式子列出来，是这样的：
其中的未知数是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 编辑 ]
页:
[1]