在微分方程计算过程中，结果中出现NaN，是怎么回事

zhuce

解微分方程运行后，出现的结果先是逐渐变大的数据，然后是变为无穷大Inf，接着就是NaN，这是怎么回事？有哪位知道啊？
7 S# N4 u. l9 t: ~: E$ P* u( E其中一段的结果如下：
& e. ]$ O( q$ k& T& Z
( Y( s$ N- U6 j& `+ ?-6.577334937091747730665057507568065230549384834272368050884038839567853*10^273,
- y; X1 u6 z, |' F* J8 P. g6 F-5.139195757419310323174537703330917750350843564381075360877238710618581*10^291，
-Inf
NaN
: O% E/ R5 \) W* D% S# INaN
NaN

  b: N4 L* w; v! [. u
是用龙格-库塔方法计算微分方程组得出的结果，程序如下：
function R = wffcz(f,g,a,b,xa,ya,N,I)
/ P4 L% x8 S. q( N2 _; w3 F%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
$ r- m* w8 X3 ~%x'=f(t,x,y) y'=g(t,x,y)
%N为迭代次数
; t* L5 `- |6 v) E/ `+ z; {! Q  X%h为步长
- e8 }" V' m# y, ^* A%ya,xa为初值
$ q6 A6 y* Q7 E' |  E* YB1=2.652e+005;A1=2.922;
B2=2.233e+005;A2=1.442;
# q/ C4 ^$ ~/ W6 _2 BR1=94.25;
C=6.897e-11;
' |/ Z2 ]4 W/ k5 Z) m* a6 ML1=21.75e-6;
2 m1 K0 f/ [' y7 p( P%f=@(t,x,y)((25079999999999999898935095034332140556544638976*t^9 - 3116000000000000083537608917399754195861504*t^8 + 164799999999999990698756673337884147712*t^7 - 4823999999999999755901658355204096*t^6 + 83320000000000003537039785984*t^5 - 770099999999999994757120*t^4 + 3330000000000000000*t^3 - 94230000000000*t^2 + 1885000000*t + 1847/50-(x-((x-B1)/A1-y)*R1-B2)/A2-(x-B1)/A1)/C);
f=@(t,x,y)((I-(x-((x-B1)/A1-y)*R1-B2)/A2-(x-B1)/A1)/C);
; T6 s5 r: x5 F4 v) z7 Ug=@(t,x,y)(((x-((x-B1)/A1-y)*R1-B2)/A2-y)*R1/L1);
h=(b-a)/N;
T=zeros(1,N+1);
X=zeros(1,N+1);
Y=zeros(1,N+1);
T=a:h:b;
6 E, L- }; w( ]4 k+ mX(1)=xa;
' o: I. F8 o+ a- ?( @Y(1)=ya;
for j=1:N
f1=vpa(feval(f,T(j),X(j),Y(j)),90);
g1=vpa(feval(g,T(j),X(j),Y(j)),90);
f2=vpa(feval(f,T(j)+h/2,X(j)+h/2*f1,Y(j)+g1/2),90);
g2=vpa(feval(g,T(j)+h/2,X(j)+h/2*f1,Y(j)+h/2*g1),90);
# X+ Z6 x, X* \8 V) Y, Sf3=vpa(feval(f,T(j)+h/2,X(j)+h/2*f2,Y(j)+h*g2/2),90);
, `7 `0 d( M9 k0 G4 lg3=vpa(feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2),90);
( w6 M7 N4 V# ^3 I. i* }, _f4=vpa(feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);
. @; [3 a4 P# h$ _g4=vpa(feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);
! p6 c, B% O6 EX(j+1)=vpa(X(j)+h*(f1+2*f2+2*f3+f4)/6,90);
' C6 f1 d* P9 Q. P6 LY(j+1)=vpa(Y(j)+h*(g1+2*g2+2*g3+g4)/6,90);
& ^, A5 y- b. V1 gR=[T' X' Y'];
end
就不知道是程序问题还是硬件问题

redplum

执行程序之前写一句“format long”,有可能是初值设的不好，另外为什么方程的系数这样大呢

athuhu

支持楼主~~~