最近在做有关ODE45的微分方程求解，在微分方程中，运用到了erf函数，虽然之前也出现过类似的问题，已经解决，但是这次丝毫找不到问题的原因所在，毫无头绪。

主程序如下：

[code]

%——————xp vp zp 分别为20个点位————————————————————————

xp=[-1.91545000000000;-1.68455000000000;-1.51545000000000;-1.28455000000000;-1.11545000000000;-0.884550000000000;-0.715450000000000;-0.484550000000000;-0.315450000000000;-0.0845500000000000;0.0845500000000000;0.315450000000000;0.484550000000000;0.715450000000000;0.884550000000000;1.11545000000000;1.28455000000000;1.51545000000000;1.68455000000000;1.91545000000000];

vp=[-2.87317500000000;-2.52682500000000;-2.27317500000000;-1.92682500000000;-1.67317500000000;-1.32682500000000;-1.07317500000000;-0.726825000000000;-0.473175000000000;-0.126825000000000;0.126825000000000;0.473175000000000;0.726825000000000;1.07317500000000;1.32682500000000;1.67317500000000;1.92682500000000;2.27317500000000;2.52682500000000;2.87317500000000];

zp=[-1.43658750000000;-1.26341250000000;-1.13658750000000;-0.963412500000000;-0.836587500000000;-0.663412500000000;-0.536587500000000;-0.363412500000000;-0.236587500000000;-0.0634125000000000;0.0634125000000000;0.236587500000000;0.363412500000000;0.536587500000000;0.663412500000000;0.836587500000000;0.963412500000000;1.13658750000000;1.26341250000000;1.43658750000000];

h0 = waitbar(0,'Formulating the Transient PDF...');

%下面进行循环，每个点位都为初始值来进行ode45的求解

for kk1=1:20

waitbar(kk1/20,h0)

for kk2=1:20

for kk3=1:20

[tt,Y] = ode45(@new_gaussisan_00way_1,[0 ,4]，[xp(kk1)^2;xp(kk1)*vp(kk2);xp(kk1)*zp(kk3);vp(kk2)^2;vp(kk2)*zp(kk3);zp(kk3)^2;xp(kk1);vp(kk2);zp(kk3)]);

end

end

end

[/code]

函数程序如下：

[code]

function dX=new_gaussisan_00way_1(t,x)

dX=zeros(9,1);

A=1;gama=0.1;beta=0.1;arf=0.1;w=1;ksai=0.05;

PI = 3.1415926;

I0 = 0.1;

z0=0;

%——上方是几个参数——————————————

%下面是几个和解有关的参数————在这些参数中，又用了了erf函数————————————————————

I101A=sqrt(2/PI)*(1/sqrt(x(4)))*(x(4)*(x(7)*x(9)+x(3))+x(2)*x(5))*exp(-x(8)^2/(2*x(4)))+(x(7)*x(5)+x(8)*x(3)+x(9)*x(2)+x(7)*x(8)*x(9))*erf(x(8)/(sqrt(2*x(4))));

I101B=sqrt(2/PI)*(1/sqrt(x(6)))*(x(6)*(x(7)*x(8)+x(2))+x(3)*x(5))*exp(-x(9)^2/(2*x(6)))+(x(7)*x(5)+x(8)*x(3)+x(9)*x(2)+x(7)*x(8)*x(9));%*erf(x(9)/(sqrt(2*x(6))));

I011A=sqrt(2/PI)*sqrt(x(4))*(x(8)*x(9)+2*x(5))*exp(-x(8)^2/(2*x(4)))+(x(9)*(x(8)^2*x(4))+2*x(8)*x(5));%*erf(x(8)/(sqrt(2*x(4))));

I011B=sqrt(2/PI)*(1/sqrt(x(6)))*(x(6)*(x(8)^2+x(4))+x(5)^2)*exp(-x(9)^2/(2*x(6)))+(x(9)*(x(8)^2*x(4))+2*x(8)*x(5));%*erf(x(9)/(sqrt(2*x(6))));

I002A=sqrt(2/PI)*(1/sqrt(x(4)))*(x(4)*(x(9)^2+x(6))+x(5)^2)*exp(-x(8)^2/(2*x(4)))+(x(8)*(x(9)^2+x(6))+2*x(5)*x(9));%*erf(x(8)/sqrt(2*x(4)));

I002B=sqrt(2/PI)*(sqrt(x(6))*(x(8)*x(9)+2*x(5)))*exp(-x(9)^2/(2*x(6)))+(x(8)*(x(9)^2+x(6))+2*x(5)*x(9));%*erf(x(9)/sqrt(2*x(6)))

I001A=sqrt(2/PI)*sqrt(x(4))*x(9)*exp(-x(8)^2/(2*x(4)))+(x(8)*x(9)+x(5));%*erf(x(8)/sqrt(2*x(4)));

I001B=sqrt(2/PI)*sqrt(x(6))*x(8)*exp(-x(9)^2/(2*x(6)))+(x(8)*x(9)+x(5));%*erf(x(9)/sqrt(2*x(6)));

%————————————————————————————————————————

dX(1)=2*x(2);

dX(2)=x(4)-2*ksai*w*x(2)-arf*w^2*x(1)-(1-arf)*w^2*x(3)+w^2*(1-arf)*z0*x(7);

dX(3)=x(5)+A*x(2)-gama*I101A-beta*I101B;

dX(4)=-4*ksai*w*x(4)-2*arf*w^2*x(2)-2*(1-arf)*w^2*x(5)+2*w^2*(1-arf)*z0*x(9)+I0*1;

dX(5)=-2*ksai*w*x(5)-arf*w^2*x(3)-(1-arf)*w^2*x(6)+w^2*(1-arf)*z0*x(9)+A*x(4)-gama*I011A-beta*I011B;

dX(6)=2*A*x(5)-2*gama*I002A-2*beta*I002B;

dX(7)=x(8);

dX(8)=-2*ksai*w*x(8)-w^2*arf*x(7)-w^2*(1-arf)*x(9)+w^2*(1-arf)*z0*I0*1;

dX(9)=A*x(8)-gama*I001A-beta*I001B;

end

[/code]

最后再运算的时候 会产生的错误如下：

[code]

错误使用 erf

输入必须为实数完全数。

出错 new_gaussisan_00way_1 (line 10)

I101A=sqrt(2/PI)*(1/sqrt(x(4)))*(x(4)*(x(7)*x(9)+x(3))+x(2)*x(5))*exp(-x(8)^2/(2*x(4)))+(x(7)*x(5)+x(8)*x(3)+x(9)*x(2)+x(7)*x(8)*x(9))*erf(x(8)/(sqrt(2*x(4))));

出错 ode45 (line 261)

f(:,2) = feval(odeFcn,t+hA(1),y+f*hB(:,1),odeArgs{:});

出错 e_rfyansuan (line 16)

[tdata,xdata]=ode45('new_gaussisan_00way_1',[0,10],[xp(kk1)^2;xp(kk1)*vp(kk2);xp(kk1)*zp(kk3);vp(kk2)^2;vp(kk2)*zp(kk3);zp(kk3)^2;xp(kk1);vp(kk2);zp(kk3)]);

[/code]

进行过一些尝试：

测试1：

如果把初始值改成0.01；0.03 这种之类的常数，就不会出现报错，

测试2：

如果初始值都采用xvp的第一个数值

[tdata,xdata]=ode45('new_gaussisan_00way_1',[0,10],[1.915^2;1.915*2.873;1.915*1.436;2.873^2;2.873*1.436;1.436^2;-1.915;-2.873;-1.436]);

也依然会使用程序的时候报错

以I101A为例，对于初始值来说erf的位置

erf(x(8)/(sqrt(2*x(4))))   就是erf(-2.873/sqrt(2*2.873^2))=erf(-1/sqrt(2))，这应该没错啊

但是依然不行。。。。

但是按道理来说，导入的初始值看起来也并没有什么问题，就是不知道哪里出错。

使得erf函数执行不下去。

希望大家能帮帮我 谢谢了