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

zhuce

解微分方程运行后，出现的结果先是逐渐变大的数据，然后是变为无穷大Inf，接着就是NaN，这是怎么回事？有哪位知道啊？
其中一段的结果如下：
% w$ p; A( \& T, _' ]$ F. q- ?0 j2 c) R
0 A0 o$ c$ Y4 ~5 G( `2 ?- Q-6.577334937091747730665057507568065230549384834272368050884038839567853*10^273,
-5.139195757419310323174537703330917750350843564381075360877238710618581*10^291，
0 w0 N# }1 v2 G# e; a2 p-Inf
/ n3 x; P, B! t  T3 Y, k- U" b* iNaN
NaN
7 c  S3 {. s$ ^/ _9 I  n$ h' P5 MNaN


是用龙格-库塔方法计算微分方程组得出的结果，程序如下：
function R = wffcz(f,g,a,b,xa,ya,N,I)
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
%x'=f(t,x,y) y'=g(t,x,y)
" e5 S4 Z8 ^, Y& w! T& }%N为迭代次数
7 ^; u) X) [/ d1 t/ b%h为步长
%ya,xa为初值
B1=2.652e+005;A1=2.922;
, ]3 {! V- ]* IB2=2.233e+005;A2=1.442;
% i- Z0 _0 L+ v* dR1=94.25;
6 }5 u# l8 h7 r5 m% ~+ nC=6.897e-11;
L1=21.75e-6;
5 d4 C; M3 n. x5 u1 w( K) N2 @% D$ r%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);
- F9 Z1 E& l) L0 z' Bf=@(t,x,y)((I-(x-((x-B1)/A1-y)*R1-B2)/A2-(x-B1)/A1)/C);
g=@(t,x,y)(((x-((x-B1)/A1-y)*R1-B2)/A2-y)*R1/L1);
5 Z( q$ Q' d+ ]2 Q9 p' {; Gh=(b-a)/N;
# l3 x4 o: k& e& n6 jT=zeros(1,N+1);
& F+ s/ c% m' Z" w5 \, x) RX=zeros(1,N+1);
5 N' V# T" _1 T6 SY=zeros(1,N+1);
T=a:h:b;
" c# _/ ^$ @/ T" W2 K4 i* o8 ~1 `X(1)=xa;
Y(1)=ya;
& H1 ~  M% _( ?8 u5 g" j* cfor j=1:N
8 K' L; y" x$ N4 K) ]9 qf1=vpa(feval(f,T(j),X(j),Y(j)),90);
- f4 P+ Y! o: f; v5 H; f8 Q3 \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);
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);
- I+ m/ G" I5 U2 v$ d- D8 ?5 A9 ag4=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);
9 p3 i- X1 n* g# a, s$ y) K1 _Y(j+1)=vpa(Y(j)+h*(g1+2*g2+2*g3+g4)/6,90);
R=[T' X' Y'];
end
就不知道是程序问题还是硬件问题

redplum

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

athuhu

支持楼主~~~