该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

在求解一微分方程时候，采用了ode45 当使用第二种初值的时候，可以完美运行求解，当使用第一种初值的时候，就会报错

报错和主程序如下， 希望大家能帮帮我

————————————————————————————————————

错误使用 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 263)

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

出错 solve_00way_1 (line 12)

[tdata,xdata]=ode45('new_gaussisan_00way_1',tspan,X0);

————————————————————————————————————

主程序如下：

clear;clc;close all

t0=cputime;

tspan=0:0.01:10;

X0=[3.6672;5.5018;2.7499;8.2541;4.1256;2.0621;-1.915;-2.873;-1.436]; 【第一种初值】

X0=[0.02,0,0,0.02,0,0.08,0.1,0.1,0.2];【第二种初值】

[tdata,xdata]=ode45('new_gaussisan_00way_1',tspan,X0);

X1=xdata(:,1);

X2=xdata(:,2);

X3=xdata(:,3);

X4=xdata(:,4);

X5=xdata(:,5);

X6=xdata(:,6);

X7=xdata(:,7);

X8=xdata(:,8);

X9=xdata(:,9);

EXX=X1;EXY=X2;EXZ=X3;EYY=X4;EYZ=X5;EZZ=X6;EX=X7;EY=X8;EZ=X9;

t1=cputime-t0;

函数程序如下：

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;

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