查看: 237|回复: 22|关注: 0

[已解决] 求解非线性方程组时计算时间太长,是不是程序有问题呢?

[复制链接]

新手

7 麦片

财富积分


050


2

主题

17

帖子

0

最佳答案
本帖最后由 李梦博1234 于 2018-12-6 11:09 编辑

非线性方程组,使用fsolve求解时一直不收敛,所以采用牛顿迭代法计算,但是计算时间太长了,程序是论坛里的,如下:
  1. function [r,m]=mulModifNewton(F,x0,h,eps)
  2. format long;
  3. if nargin==3
  4. eps=1.0e-4;
  5. end

  6. n = length(x0);
  7. fx = subs(F,findsym(F),x0);
  8. J = zeros(n,n);
  9. for i=1:n
  10. x1 = x0;
  11. x1(i) = x1(i)+h(i);
  12. J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);
  13. end
  14. r=transpose(x0)-inv(J)*fx;
  15. m=1;
  16. tol=1;
  17. while tol>eps
  18. xs=r;
  19. fx = subs(F,findsym(F),xs);
  20. J = zeros(n,n);
  21. h = fx;
  22. for i=1:n
  23. x1 = xs;
  24. x1(i) = x1(i)+h(i);
  25. J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);
  26. end
  27. r=xs-inv(J)*fx; %核心迭代公式
  28. tol=norm(r-xs);
  29. m=m+1;
  30. if(m>100000) %迭代步数控制
  31. disp('迭代步数太多,可能不收敛!');
  32. return;
  33. end
  34. end
  35. format short;
复制代码

方程及运行程序为:
  1. syms Ga T1 T3 T5 lamd w Qh
  2. z = [(1-Ga/500)^(-2)*log((Ga*(T3-306))/(Ga*(T3-303)-1500))+(1/Ga-1/500)^(-1)*(1/Ga-(T3-303)/(Ga*(T3-303)-1500))-(Ga/400-1)^(-2)*log((Ga*(293-T1))/(Ga*(295.5-T1)+1000))+(1/400-1/Ga)^(-1)*lamd*(T1/(Ga*T1-1000)-(T1+T3-T5)/(Ga*(T1+T3-T5)-1500))+w*(T3-T5);
  3. -(1/400-1/Ga)^(-1)*(-Ga/(Ga*(295.5-T1)+1000)+1/(293-T1))+(Ga*(T3-T5)-1500)/(T5-T1)^2+lamd*Ga*(1/(Ga*T1-1000)-1/(Ga*(T1+T3-T5)-1500))+w*10;
  4. (1/Ga-1/500)^(-1)*(1/(T3-306)-Ga/(Ga*(T3-303)-1500))+Ga/(T5-T1)-lamd*(Ga/(Ga*(T1+T3-T5)-1500)-1/T3)+w*Ga;
  5. (Ga*(T1-T3)+1500)/(T5-T1)^2+lamd*(-1/T5+Ga/(Ga*(T3+T1-T5)-1500))-w*(Ga+10);
  6. Ga*(T3-T5)-1500-10*(T5-T1);
  7. log((Ga*T1-1000)/T5)-log((Ga*(T1+T3-T5)-1500)/T3)];
  8. [r,m] = mulModifNewton(z,([1 2 3 4 5 6]),[0.1 0.1 0.1 0.1 0.1 0.1])
复制代码
六个方程组,使用fsolve一直不收敛,exitflag=0。我使用MATLAB2017a,电脑四核8G运存。

还请各位大神帮忙看看,万分感谢!

QQ截图20181206105457.jpg

新手

7 麦片

财富积分


050


2

主题

17

帖子

0

最佳答案
 楼主| 发表于 7 天前 | 显示全部楼层
  1. function [r,m]=mulModifNewton(F,x0,h,eps)
  2. format long;
  3. if nargin==3
  4. eps=1.0e-4;
  5. end

  6. n = length(x0);
  7. fx = subs(F,findsym(F),x0);
  8. J = zeros(n,n);
  9. for i=1:n
  10. x1 = x0;
  11. x1(i) = x1(i)+h(i);
  12. J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);
  13. end
  14. r=transpose(x0)-inv(J)*fx;
  15. m=1;
  16. tol=1;
  17. while tol>eps
  18. xs=r;
  19. fx = subs(F,findsym(F),xs);
  20. J = zeros(n,n);
  21. h = fx;
  22. for i=1:n
  23. x1 = xs;
  24. x1(i) = x1(i)+h(i);
  25. J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);
  26. end
  27. r=xs-inv(J)*fx; %核心迭代公式
  28. tol=norm(r-xs);
  29. m=m+1;
  30. if(m>100000) %迭代步数控制
  31. disp('迭代步数太多,可能不收敛!');
  32. return;
  33. end
  34. end
  35. format short;
复制代码

新手

7 麦片

财富积分


050


2

主题

17

帖子

0

最佳答案
 楼主| 发表于 7 天前 | 显示全部楼层

这个是程序代码,上面出错了,没有贴出来,麻烦大神看看呀,谢谢

论坛优秀回答者

中级

851 麦片

财富积分


5001500


0

主题

2172

帖子

185

最佳答案
  • 关注者: 117
发表于 7 天前 | 显示全部楼层
参考下1stOpt的计算结果:
ga: 3.47632463709284E-8
t3: -503.944746269951
t1: -185.963163012407
lamd: -2.3086997261903E-8
t5: -335.963163596366
w: 0.00666666663366633

新手

7 麦片

财富积分


050


2

主题

17

帖子

0

最佳答案
 楼主| 发表于 7 天前 | 显示全部楼层
shihe 发表于 2018-12-6 11:16
参考下1stOpt的计算结果:
ga: 3.47632463709284E-8
t3: -503.944746269951

非常感谢您,但是结果不太对

论坛优秀回答者

中级

851 麦片

财富积分


5001500


0

主题

2172

帖子

185

最佳答案
  • 关注者: 117
发表于 7 天前 | 显示全部楼层
数值结果的话哪儿不对?
ga: 3.74740299151007E-10
t3: -503.944752866633
t1: -185.963168565183
lamd: -2.39095681275822E-10
t5: -335.963168571478
w: 0.00666666666631295

新手

7 麦片

财富积分


050


2

主题

17

帖子

0

最佳答案
 楼主| 发表于 7 天前 | 显示全部楼层
shihe 发表于 2018-12-6 14:00
数值结果的话哪儿不对?
ga: 3.74740299151007E-10
t3: -503.944752866633

ga是22,温度都是正的,270左右,400左右,300左右,您这个算出来怎么都是负数呢?

论坛优秀回答者

中级

851 麦片

财富积分


5001500


0

主题

2172

帖子

185

最佳答案
  • 关注者: 117
发表于 7 天前 | 显示全部楼层
从数学角度6#的结果没任何问题,代入验证下就知道了。至于你提到的参数范围要求应该提前说明。
6个方程6个参数,要满足你范围要求的精确解估计不存在。

新手

7 麦片

财富积分


050


2

主题

17

帖子

0

最佳答案
 楼主| 发表于 7 天前 | 显示全部楼层
shihe 发表于 2018-12-6 16:55
从数学角度6#的结果没任何问题,代入验证下就知道了。至于你提到的参数范围要求应该提前说明。
6个方程6个 ...

非常感谢您的回答,需要达到我的那组范围解,应该需要几组方程呢?

论坛优秀回答者

中级

851 麦片

财富积分


5001500


0

主题

2172

帖子

185

最佳答案
  • 关注者: 117
发表于 7 天前 | 显示全部楼层
你最好给出每个参数具体的范围限制
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

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