声振论坛

 找回密码
 我要加入

QQ登录

只需一步,快速开始

查看: 3842|回复: 8

[编程技巧] 如何解决程序循环次数增加后出现Inf或NAN

[复制链接]
发表于 2007-7-21 10:39 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

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

x
:@( 在matlab矩阵计算
X=A/B
前300个循环一切正常
300以后突然发散成Inf或NAN
可能是甚么原因啊?
麻烦高手指点

[ 本帖最后由 eight 于 2007-7-23 12:50 编辑 ]
回复
分享到:

使用道具 举报

发表于 2007-7-21 10:47 | 显示全部楼层
能不能把程序写多点啊!A B 是什么啊?
 楼主| 发表于 2007-7-21 19:29 | 显示全部楼层
本帖最后由 牛小贱 于 2014-6-29 13:14 编辑

AX=B
A是系数矩阵,B是列矩阵
附程序如下,望各位高手不吝赐教,小女子暂且先谢过:
  1. clear
  2. yy=1e-6;zz=1e-6;tt=1e-3;yita=4.6857e-7; R=3.3333e-5; seta=4.027e+6; B=0.5;rou=2450;
  3. N=1;
  4. A=5e-4;
  5. m=200;n=500;
  6. C=zeros(2*m+1,3*n); E=zeros(2*m+1,3*n);
  7. F=zeros(2*m+1,3*n); H=zeros(2*m+1,3*n);
  8. D=939140;G=-935140;
  9. X=zeros(2*m+1,3*n)+eps;
  10. QX=zeros(2*m+1,3*n)+eps;
  11. U=ones(2*m+1,3*n)*eps;
  12. QU=ones(2*m+1,3*n)*eps;

  13. vy=zeros(2*m+1,3*n);
  14. vz=zeros(2*m+1,3*n);
  15. C=-[(N*vy)/(2*yy)]-yita/yy^2; E=[(N*vy)/(yy^2)]-yita/yy^2;
  16. F=[(N*(vz+R))/(2*zz)]+yita/zz^2;
  17. H=-[(N*(vz+R))/(2*zz)]+yita/zz^2;


  18. A1=zeros(2*m+1);
  19. for j=2:n-1
  20. i=m;
  21. while i*yy>800*(j*zz)^2
  22. i=i-1;
  23. end

  24.                                                 
  25. A1(1+m-i,1+m-i)=1;
  26. A1(1+m+i,1+m+i)=1;
  27. for k=2+m-i:m+i
  28. A1(k,k-1)=C(k,k-1);A1(k,k)=D;A1(k,k+1)=E(k,k+1);
  29. end

  30. AA=zeros(2*i+1);
  31. u=1;v=1;
  32. for h=1+m-i:1+m+i
  33. for g=1+m-i:1+m+i
  34. AA(u,v)=A1(h,g);
  35. v=v+1;
  36. end
  37. u=u+1;
  38. v=1;
  39. end

  40. Y=zeros(2*m+1,1);
  41. for f=2+m-i:m+i
  42. Y(f,1)=F(f,j-1)*X(f,j)+G*X(f,j)+H(f,j+1)*X(f,j+1)-seta*(B^2)*X(f,j)/rou;
  43. end
  44. Y(1+m-i,1)=2*eps/yy^2;
  45. Y(1+m+i,1)=2*eps/yy^2;

  46. Y1=zeros(2*i+1,1);
  47. u=1;
  48. for h=1+m-i:1+m+i
  49. Y1(u,1)=Y(h,1);
  50. u=u+1;
  51. end

  52. X1=zeros(2*i+1,1);
  53. X1=(AA+eps)\Y1;

  54. Xj=zeros(2*m+1,1);
  55. u=1;
  56. for h=1+m-i:1+m+i
  57. Xj(h,1)=X1(u,1);
  58. u=u+1;
  59. end

  60. X(:,j)=Xj;
  61. end                                          
  62. X
复制代码


发表于 2007-7-21 22:12 | 显示全部楼层
试试pinv.

评分

1

查看全部评分

 楼主| 发表于 2007-7-22 09:49 | 显示全部楼层
改成
X1=pinv((AA+eps)\Y1);
运行出现:

??? Index exceeds matrix dimensions
Error in   ==> Xj(h,1)=X1(u,1);

pinv试了 不好用也:@(

[ 本帖最后由 ChaChing 于 2010-4-4 20:57 编辑 ]
发表于 2007-7-23 12:49 | 显示全部楼层
请看看置顶贴,然后找到“常见出错问题整理”的帖子,建议先阅读一下

[ 本帖最后由 ChaChing 于 2010-4-4 20:59 编辑 ]

评分

1

查看全部评分

发表于 2011-8-7 15:27 | 显示全部楼层
回复 1 # huangxuemei513 的帖子

请问,您的问题得到解决了么?我现在也遇到了同样的问题,计算过程总是应为NaN终止,比如100次的迭代,一般到10次以下就终止了,望指点
发表于 2014-6-19 21:44 | 显示全部楼层
huangxuemei513 发表于 2007-7-21 19:29
AX=B
A是系数矩阵,B是列矩阵
附程序如下,望各位高手不吝赐教,小女子暂且先谢过:

你好,你的问题解决了吗?我在分解循环的时候也遇到矩阵的最后一行是INF,望指点,谢谢!
发表于 2014-6-20 09:09 | 显示全部楼层
liyukui89 发表于 2014-6-19 21:44
你好,你的问题解决了吗?我在分解循环的时候也遇到矩阵的最后一行是INF,望指点,谢谢!

参考4楼方法,试试pinv
您需要登录后才可以回帖 登录 | 我要加入

本版积分规则

QQ|小黑屋|Archiver|手机版|联系我们|声振论坛

GMT+8, 2024-12-23 19:56 , Processed in 0.084328 second(s), 19 queries , Gzip On.

Powered by Discuz! X3.4

Copyright © 2001-2021, Tencent Cloud.

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