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

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

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

[复制链接]

该用户从未签到

尚未签到

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

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

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

×
解微分方程运行后,出现的结果先是逐渐变大的数据,然后是变为无穷大Inf,接着就是NaN,这是怎么回事?有哪位知道啊?
# B) i; r% k( }& r% E其中一段的结果如下:) g9 Q! @0 S, R) r' p7 y
' c' ?& E" ^$ }$ T" C
-6.577334937091747730665057507568065230549384834272368050884038839567853*10^273,
, F0 [, a* Q7 p4 {7 V; m-5.139195757419310323174537703330917750350843564381075360877238710618581*10^291,/ K: J& E9 a& k3 V9 r4 V" w
-Inf
2 r* L, O8 ]  r) J" |0 |% A' eNaN
7 r/ d( K% _8 s& QNaN
  ^8 E/ Y+ Y# O% q( B2 _: H8 lNaN6 C& \; u$ {; i9 m. M; u7 R
5 L8 o4 N+ a! n3 F% m
- M3 U. b% x* _$ B, V
是用龙格-库塔方法计算微分方程组得出的结果,程序如下:
& W" F6 \+ J$ l7 o8 v8 xfunction R = wffcz(f,g,a,b,xa,ya,N,I)7 T" h/ {: C1 l; b
%UNTITLED2 Summary of this function goes here- h4 T/ Y, w0 C: Z2 Z3 ~/ u
% Detailed explanation goes here0 D/ n1 K6 {. D2 W0 I0 D" m
%x'=f(t,x,y) y'=g(t,x,y)4 N- M' G7 T0 \
%N为迭代次数7 e  r8 z% u" d' C3 H% c- G
%h为步长
; a0 m5 K/ q5 P& x" w1 C: ?2 N%ya,xa为初值8 y6 J6 O+ k+ h: f& i% P$ A- G8 E6 D
B1=2.652e+005;A1=2.922;
0 h$ X6 {" M3 B* B/ lB2=2.233e+005;A2=1.442;
* M( b3 ^7 F$ t( \R1=94.25;
& R% a; P7 T3 b- ~9 G. rC=6.897e-11;7 y$ T; b+ t- b) C! E. ]6 Q
L1=21.75e-6;
) ~; z; K0 v2 A: M%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);
) l! U4 h" E  F* R! j  k% Mf=@(t,x,y)((I-(x-((x-B1)/A1-y)*R1-B2)/A2-(x-B1)/A1)/C);7 G( L1 ^9 J8 O
g=@(t,x,y)(((x-((x-B1)/A1-y)*R1-B2)/A2-y)*R1/L1);% t1 N1 G% Q* \! `$ d# ~6 t
h=(b-a)/N;
4 k0 T0 w8 ]# e/ M0 @6 E- fT=zeros(1,N+1);
0 P8 \* j( g5 r: s# FX=zeros(1,N+1);9 F7 u8 q( {& L2 B) k9 v1 N
Y=zeros(1,N+1);5 V" A5 T: i$ F" @3 x
T=a:h:b;
& s2 j0 c" R- A1 V) F' w# O: FX(1)=xa;
0 @8 a8 K+ w" Q  X, n: @Y(1)=ya;
' f/ _+ H8 W$ t4 ofor j=1:N
; K  f+ d4 w, K4 W! u; J, n+ |f1=vpa(feval(f,T(j),X(j),Y(j)),90);
) \; u3 B! n, f$ Jg1=vpa(feval(g,T(j),X(j),Y(j)),90);' m& A9 A3 M1 v  D4 l/ k8 q
f2=vpa(feval(f,T(j)+h/2,X(j)+h/2*f1,Y(j)+g1/2),90);, ~/ p/ r/ M+ M) N8 \9 C/ s# k: N
g2=vpa(feval(g,T(j)+h/2,X(j)+h/2*f1,Y(j)+h/2*g1),90);
$ ^! h; x# G  \f3=vpa(feval(f,T(j)+h/2,X(j)+h/2*f2,Y(j)+h*g2/2),90);4 z1 k% e6 P9 z/ p0 v7 m
g3=vpa(feval(g,T(j)+h/2,X(j)+h/2*f2,Y(j)+h/2*g2),90);5 G- M6 `- p* O5 @3 F
f4=vpa(feval(f,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);
9 L3 b# T, |  s' g. n- j: tg4=vpa(feval(g,T(j)+h,X(j)+h*f3,Y(j)+h*g3),90);6 w0 @1 W  @* ~9 H  o
X(j+1)=vpa(X(j)+h*(f1+2*f2+2*f3+f4)/6,90);. Y" h+ ]* W# u
Y(j+1)=vpa(Y(j)+h*(g1+2*g2+2*g3+g4)/6,90);: [0 @1 O  @1 T( A+ w
R=[T' X' Y'];
/ o& i. m5 N& X( j' K- jend
* p4 S$ y: _, ?  L0 L/ A就不知道是程序问题还是硬件问题
楼主热帖
帖文化:【文明发帖 和谐互动】 社区精神:【创新、交流、互助、共享】
  • 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 18:38

    Powered by Discuz! X3.5 Licensed

    © 2001-2025 Discuz! Team.

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