最近在做有关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函数执行不下去。
希望大家能帮帮我 谢谢了