设为首页收藏本站|繁體中文 快速切换版块

 找回密码
 立即加入
搜索
查看: 963|回复: 2

[讨论] 在微分方程计算过程中,结果中出现NaN,是怎么回事

[复制链接]

该用户从未签到

尚未签到

发表于 2013-2-28 17:14:03 | 显示全部楼层 |阅读模式

马上加入,结交更多好友,共享更多资料,让你轻松玩转电力研学社区!

您需要 登录 才可以下载或查看,没有账号?立即加入

×
解微分方程运行后,出现的结果先是逐渐变大的数据,然后是变为无穷大Inf,接着就是NaN,这是怎么回事?有哪位知道啊?
- Y7 W: C( N, q" ]) ~" I6 m7 c其中一段的结果如下:/ j7 D, e- K6 {5 a; }2 G0 E6 a
3 w5 S7 U3 G: g
-6.577334937091747730665057507568065230549384834272368050884038839567853*10^273,6 T3 ]% H) U% _9 h' M
-5.139195757419310323174537703330917750350843564381075360877238710618581*10^291,
# H2 W5 h- ^/ ?* t8 b-Inf5 H9 Q" Q7 u: z. F
NaN/ I2 O2 g$ f5 ^* |- P- b! B& F
NaN- J# l/ C9 i0 Y+ l$ k
NaN, j. X3 m+ _. _/ ^

. L9 e- u, ^% `$ C" w$ t- S2 z! _; H! R2 i3 V
是用龙格-库塔方法计算微分方程组得出的结果,程序如下:! z% c' G. [% z; Q/ r
function R = wffcz(f,g,a,b,xa,ya,N,I)! {& @+ i7 a& Y- r! I4 ^" D7 }5 h" e) }
%UNTITLED2 Summary of this function goes here2 u$ }3 l3 v" i9 z  e$ G/ J
% Detailed explanation goes here5 U  ^7 V( Z7 d* K/ Q9 e) c5 @
%x'=f(t,x,y) y'=g(t,x,y)2 }6 R1 B' g6 b
%N为迭代次数
8 h) J* @0 a8 a8 q3 ~$ u%h为步长
0 a  K  a7 R7 y  b%ya,xa为初值
$ J3 R: d& S& d9 h. I1 xB1=2.652e+005;A1=2.922;
! H0 O$ a- w- g$ H$ c; f0 a6 dB2=2.233e+005;A2=1.442;0 ], s/ w8 `; v$ B: F2 {% Q
R1=94.25;
6 s! d- p2 L) n0 ~# @! M! `; w& bC=6.897e-11;$ R# ~/ x; Y. Y/ V4 I
L1=21.75e-6;' B7 I. g3 O! k* H
%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);0 A( S! d  _' r$ m9 Q$ c( u  B1 O
f=@(t,x,y)((I-(x-((x-B1)/A1-y)*R1-B2)/A2-(x-B1)/A1)/C);
  c+ {2 }; Q9 Q$ m* M) {+ ~g=@(t,x,y)(((x-((x-B1)/A1-y)*R1-B2)/A2-y)*R1/L1);
% X! ]2 v" i, D- N6 oh=(b-a)/N;$ k1 \6 X( Z9 a; u! Z
T=zeros(1,N+1);
) Z1 g/ f6 @& v4 I; Y0 ]4 LX=zeros(1,N+1);
- @* L( v5 X1 l" D& i/ SY=zeros(1,N+1);; Y6 v( ~5 K# \+ I, g
T=a:h:b;- Q- L  q) s9 u8 P
X(1)=xa;% n' ]  g  q1 F0 Y+ i% E+ R5 x8 [
Y(1)=ya;& E: W% a3 V) S9 y, q
for j=1:N/ v7 _# @. a; y/ [) b/ f& [
f1=vpa(feval(f,T(j),X(j),Y(j)),90);
7 g& n- J, H. J+ [4 Eg1=vpa(feval(g,T(j),X(j),Y(j)),90);+ ?" u5 m7 R# @, i
f2=vpa(feval(f,T(j)+h/2,X(j)+h/2*f1,Y(j)+g1/2),90);- O: s$ z7 b, A. @' P
g2=vpa(feval(g,T(j)+h/2,X(j)+h/2*f1,Y(j)+h/2*g1),90);
' Q/ n: f/ }* U3 _, @; Df3=vpa(feval(f,T(j)+h/2,X(j)+h/2*f2,Y(j)+h*g2/2),90);: w6 r: `& T- B
g3=vpa(feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2),90);
5 H, n2 e9 |! i/ O1 }" y4 [7 Bf4=vpa(feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);! S" `1 v. n( _- A  \! L+ s+ }- j
g4=vpa(feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);
) x3 O( m3 l3 X- Q+ b5 _, rX(j+1)=vpa(X(j)+h*(f1+2*f2+2*f3+f4)/6,90);
* ?, G& p0 \4 _Y(j+1)=vpa(Y(j)+h*(g1+2*g2+2*g3+g4)/6,90);0 Y9 b/ K$ P0 d: i1 ^* F$ d
R=[T' X' Y'];
) m! J  o7 @$ ~# Y2 iend
* a# A, l* V2 O就不知道是程序问题还是硬件问题
楼主热帖
帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
  • TA的每日心情
    愤怒
    2021-6-12 00:00
  • 签到天数: 1657 天

    连续签到: 28 天

    [LV.Master]伴坛终老

    累计签到:2776 天
    连续签到:15 天
    发表于 2013-3-1 00:36:11 | 显示全部楼层
    执行程序之前写一句“format long”,有可能是初值设的不好,另外为什么方程的系数这样大呢
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】

    该用户从未签到

    尚未签到

    发表于 2013-3-1 08:23:30 | 显示全部楼层
    支持楼主~~~
    "真诚赞赏,手留余香"
    还没有人打赏,支持一下
    帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
    您需要登录后才可以回帖 登录 | 立即加入

    本版积分规则

    招聘斑竹

    小黑屋|手机版|APP下载(beta)|Archiver|电力研学网 ( 赣ICP备12000811号-1|赣公网安备36040302000210号 )|网站地图

    GMT+8, 2025-4-10 14:10

    Powered by Discuz! X3.5 Licensed

    © 2001-2025 Discuz! Team.

    快速回复 返回顶部 返回列表