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

zhuce

解微分方程运行后，出现的结果先是逐渐变大的数据，然后是变为无穷大Inf，接着就是NaN，这是怎么回事？有哪位知道啊？
3 Q$ u# k1 a# P$ M# i其中一段的结果如下：
  l2 I* t5 M0 O5 X5 D
-6.577334937091747730665057507568065230549384834272368050884038839567853*10^273,
-5.139195757419310323174537703330917750350843564381075360877238710618581*10^291，
-Inf
7 r! q4 O1 R$ r$ ~8 w) x! YNaN
NaN
NaN


是用龙格-库塔方法计算微分方程组得出的结果，程序如下：
  o8 s7 z" F: K  |" Qfunction R = wffcz(f,g,a,b,xa,ya,N,I)
% g+ ~, h/ o, S& }4 C! r%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
$ k# }1 N/ |5 E3 D%x'=f(t,x,y) y'=g(t,x,y)
%N为迭代次数
%h为步长
%ya,xa为初值
B1=2.652e+005;A1=2.922;
4 E- O4 f9 C* f/ y  b1 _B2=2.233e+005;A2=1.442;
R1=94.25;
) d% x- j/ U1 ^# ?  BC=6.897e-11;
" h$ l+ K1 G; x* iL1=21.75e-6;
$ K- P# D& }; T%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);
% r% x" ~+ w8 J% ]: j2 Z0 C4 |g=@(t,x,y)(((x-((x-B1)/A1-y)*R1-B2)/A2-y)*R1/L1);
h=(b-a)/N;
* F) ]" y% q, L7 K1 N) x% pT=zeros(1,N+1);
* `( y! G. M0 V7 x6 |; f0 aX=zeros(1,N+1);
9 ~1 G+ F. j" q" |! h; IY=zeros(1,N+1);
8 @3 q$ p. U% U7 S7 ~T=a:h:b;
. ]* ^" N- @7 V6 PX(1)=xa;
: ^5 F6 o  d6 R% w5 _7 OY(1)=ya;
# {' o7 s$ D* C  |' ofor 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);
' _5 U* A: w5 ?/ K' If2=vpa(feval(f,T(j)+h/2,X(j)+h/2*f1,Y(j)+g1/2),90);
0 X. O( g' ?  ^  v6 lg2=vpa(feval(g,T(j)+h/2,X(j)+h/2*f1,Y(j)+h/2*g1),90);
f3=vpa(feval(f,T(j)+h/2,X(j)+h/2*f2,Y(j)+h*g2/2),90);
g3=vpa(feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2),90);
f4=vpa(feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);
g4=vpa(feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);
X(j+1)=vpa(X(j)+h*(f1+2*f2+2*f3+f4)/6,90);
( J, T$ T6 A& hY(j+1)=vpa(Y(j)+h*(g1+2*g2+2*g3+g4)/6,90);
* _& ], m* I3 o$ p' kR=[T' X' Y'];
! P* Y; @* I# N' X3 f; Z( yend
2 @2 R9 ~# Y( c- Q( H6 X/ `% ~就不知道是程序问题还是硬件问题

redplum

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

athuhu

支持楼主~~~