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

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

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

[复制链接]

该用户从未签到

尚未签到

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

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

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

×
解微分方程运行后,出现的结果先是逐渐变大的数据,然后是变为无穷大Inf,接着就是NaN,这是怎么回事?有哪位知道啊?
7 ]: i; H5 _1 |) w) T: f其中一段的结果如下:* ~- Y+ ?. l3 R0 F  A: J" \. f- y
/ {$ Y; v' K% _
-6.577334937091747730665057507568065230549384834272368050884038839567853*10^273,
4 ?  Y0 y) o3 h5 A! _-5.139195757419310323174537703330917750350843564381075360877238710618581*10^291,5 n; f! f( L3 J. k
-Inf
4 s! M& A# Z; j$ K: mNaN9 v, q% ^: ?. F4 e% p$ O% H. U& l
NaN
6 _2 [( |. L  [" i2 F# M# p# tNaN
3 T2 M5 X& q' m
! X: ]) Q( @  F' |, j% L7 m5 i0 g, e, {" S" u# m8 }; w
是用龙格-库塔方法计算微分方程组得出的结果,程序如下:
4 |; H9 ~+ V2 G) f5 yfunction R = wffcz(f,g,a,b,xa,ya,N,I)  U+ J2 T2 D3 D% Z2 V
%UNTITLED2 Summary of this function goes here
) v" \5 {8 W, m% Detailed explanation goes here3 z% r/ f) N* C6 ?
%x'=f(t,x,y) y'=g(t,x,y)
5 {! ]; q/ z  z) j% ]- M%N为迭代次数
3 I% P# S7 i6 F! f4 j%h为步长. V; Q  Z, k0 K+ L4 {2 q7 n; D5 g
%ya,xa为初值- Q& Y% i" b& `: [
B1=2.652e+005;A1=2.922;9 @. _+ h$ ]' Q5 a( ~( Z
B2=2.233e+005;A2=1.442;
$ z9 ^5 E  z/ O" r. GR1=94.25;% [, G) F# X  F- F; J/ ^0 M3 v6 q0 \
C=6.897e-11;" [* K1 f! k( {7 W, o
L1=21.75e-6;
( x5 x2 P8 {1 g, m; Y' `3 A%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);
/ U1 J' K- }4 B3 N" m* e- \" T# I% Qf=@(t,x,y)((I-(x-((x-B1)/A1-y)*R1-B2)/A2-(x-B1)/A1)/C);
+ _* q! ]) M- D# Jg=@(t,x,y)(((x-((x-B1)/A1-y)*R1-B2)/A2-y)*R1/L1);
- ?+ C5 [" M9 [  Dh=(b-a)/N;
3 `" P( p) B* M/ c" T: MT=zeros(1,N+1);3 I; B$ g6 T2 a0 x) ~1 Z
X=zeros(1,N+1);4 F! l  [* I6 j( C! V. E
Y=zeros(1,N+1);8 V, }  k0 n" s1 \
T=a:h:b;
/ i; R2 X$ ?) y: a) |* b$ PX(1)=xa;, P( M' [* C0 A5 f" F+ E- u! ?1 _
Y(1)=ya;
- `' @8 _" C( a) L% S, G  G9 q) wfor j=1:N2 S% I6 D  j- M  B8 u2 R
f1=vpa(feval(f,T(j),X(j),Y(j)),90);
) J7 Z- K6 C! D- Qg1=vpa(feval(g,T(j),X(j),Y(j)),90);
+ H3 k7 V7 f* u3 ?6 ?! Ef2=vpa(feval(f,T(j)+h/2,X(j)+h/2*f1,Y(j)+g1/2),90);
4 n& s# ~9 [  `" s7 u0 O8 Ug2=vpa(feval(g,T(j)+h/2,X(j)+h/2*f1,Y(j)+h/2*g1),90);
( n0 E8 W3 b7 `  M4 U% @& h5 P6 h/ S. Df3=vpa(feval(f,T(j)+h/2,X(j)+h/2*f2,Y(j)+h*g2/2),90);0 T* o, c) @! M6 |) q/ L. ]7 p
g3=vpa(feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2),90);
8 y7 V3 C! `( A: l1 x1 t- Y1 yf4=vpa(feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);; k( n9 G; O# Y  E' m
g4=vpa(feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);# K; |" J  l8 @
X(j+1)=vpa(X(j)+h*(f1+2*f2+2*f3+f4)/6,90);0 m( k& e( N- E
Y(j+1)=vpa(Y(j)+h*(g1+2*g2+2*g3+g4)/6,90);7 a6 V) _+ D' ]
R=[T' X' Y'];
9 `: x3 _7 R' C5 r  c) Send9 G3 r( \! r# [6 m+ q: S; E- U
就不知道是程序问题还是硬件问题
楼主热帖
帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
  • TA的每日心情
    愤怒
    2021-6-12 00:00
  • 签到天数: 1657 天

    连续签到: 28 天

    [LV.Master]伴坛终老

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

    该用户从未签到

    尚未签到

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

    本版积分规则

    招聘斑竹

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

    GMT+8, 2026-3-16 13:40

    Powered by Discuz! X3.5 Licensed

    © 2001-2025 Discuz! Team.

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