[已解决] 牛顿迭代求面积,但是结果出来的永远是负数和虚数。

[复制链接]
wikys 发表于 2021-11-20 07:10:57
图片是个物理算式,不过不重要,计算导热管的面积。所有数据都给了,只有d是未知量。

我的代码运行出来的最后一直是负数和虚数,(面积起码得是个正数,所以我知道哪里一定写错了)不太懂问题出在哪儿。求大佬带带我:'(

我的思路是:转写图片上的方程,再求导,使 f = 0,然后再用牛顿法的式子: x - f/f‘。

数字比较多。。。辛苦各位了!!感恩的心

% di = 1.219;
% Rfi = 1.76 * 10^(-4);
% Rfo = 1.76 * 10^(-4);
% hs = 356;
% ht = 356;
% kw = 60;
%
deltaT = 29.6;
Q = 801368;
Aex = 64.15;


x = 0.007;
U1 = x/(1.219 * 356);
U2 = (x * 1.76 * 10^(-4))/1.219;
U3 = (x/1.219);
U4 = x * log(U3)/(2*60);
U5 = 1/(356);
U6 = 1.76 * 10^(-4);
Uf = 1./(U1+U2+U6+U4+U5);
% a = 1.219;
% b = 356;
% c = 1.76 * 10^(-4);
% d = 1.76 * 10^(-4);
% e = 356;
% T = 29.6;
% Q = 801368;
U7 = 801368/(29.6*1.219*356);
U8 = 801368/(120*29.6);
U9 = (801368*1.76 * 10^(-4))/(29.6*1.219);
U10 = (801368*log(x/1.219))/(120*29.6)
f = @(x) (Q /(Uf * deltaT)) - Aex;
fprim = @(x) U7+U8+U9+U10;
% x0 = 0.007;
% tol = 10^(-8);
% max = 30;
% x(1) = d0;
x = 0.007;
x_old = 100;
x_true = 10^(-8);
iter = 0;
while abs(x_old-x) > 10^-3 && x ~= 0
    x_old = x;
    x = x - ((Q /(Uf * deltaT)) - Aex)/(U7+U8+U9+U10);
    iter = iter + 1;
    fprintf('Iteration %d: x=%.20f, err=%.20f\n',iter, x, x_true-x);
    pause;
end


截屏2021-11-19 23.51.37.png

最佳答案


TouAkira 发表于 2021-11-20 07:58:41
仔细看循环,里面每轮迭代的更新,根本就没有使用到 f / f‘ 在当前x值时的取值,也就是说,在最关键的地方刻舟求剑了。
你用的  x = x - ((Q /(Uf * deltaT)) - Aex)/(U7+U8+U9+U10) 标粗的部分都是你前面算出来的具体数值,而不是能够根据x值变化而更新数值的函数
  1. clear;
  2. syms x real;
  3. deltaT = 29.6;
  4. Q = 801368;
  5. Aex = 64.15;
  6. U1 = x/(1.219 * 356);
  7. U2 = (x * 1.76 * 10^(-4))/1.219;
  8. U3 = (x/1.219);
  9. U4 = x * log(U3)/(2*60);
  10. U5 = 1/(356);
  11. U6 = 1.76 * 10^(-4);
  12. Uf = 1./(U1+U2+U6+U4+U5);
  13. U7 = 801368/(29.6*1.219*356);
  14. U8 = 801368/(120*29.6);
  15. U9 = (801368*1.76 * 10^(-4))/(29.6*1.219);
  16. U10 = (801368*log(x/1.219))/(120*29.6);
  17. f = (Q /(Uf * deltaT)) - Aex;
  18. fprim = diff( f, x, 1 );
  19. f = matlabFunction( f, 'vars', x ); % 这样才是函数!
  20. fprim = matlabFunction( fprim, 'vars', x );
  21. clear x;
  22. x = 0.007;
  23. iter = 0;
  24. while abs( f( x ) ) > 10^-6 && x ~= 0
  25.     x = x - f( x )/ fprim( x ); % 使用函数,根据前一步x的数值去计算对应的f和导数fprim,然后更新
  26.     iter = iter + 1;
  27.     fprintf('Iteration %d: x=%.20f, err=%.20f\n',iter, x, abs( f( x ) ) );
  28.     pause(.1);
  29. end
复制代码
回复此楼

2 条回复


TouAkira 发表于 2021-11-20 07:58:41
仔细看循环,里面每轮迭代的更新,根本就没有使用到 f / f‘ 在当前x值时的取值,也就是说,在最关键的地方刻舟求剑了。
你用的  x = x - ((Q /(Uf * deltaT)) - Aex)/(U7+U8+U9+U10) 标粗的部分都是你前面算出来的具体数值,而不是能够根据x值变化而更新数值的函数
  1. clear;
  2. syms x real;
  3. deltaT = 29.6;
  4. Q = 801368;
  5. Aex = 64.15;
  6. U1 = x/(1.219 * 356);
  7. U2 = (x * 1.76 * 10^(-4))/1.219;
  8. U3 = (x/1.219);
  9. U4 = x * log(U3)/(2*60);
  10. U5 = 1/(356);
  11. U6 = 1.76 * 10^(-4);
  12. Uf = 1./(U1+U2+U6+U4+U5);
  13. U7 = 801368/(29.6*1.219*356);
  14. U8 = 801368/(120*29.6);
  15. U9 = (801368*1.76 * 10^(-4))/(29.6*1.219);
  16. U10 = (801368*log(x/1.219))/(120*29.6);
  17. f = (Q /(Uf * deltaT)) - Aex;
  18. fprim = diff( f, x, 1 );
  19. f = matlabFunction( f, 'vars', x ); % 这样才是函数!
  20. fprim = matlabFunction( fprim, 'vars', x );
  21. clear x;
  22. x = 0.007;
  23. iter = 0;
  24. while abs( f( x ) ) > 10^-6 && x ~= 0
  25.     x = x - f( x )/ fprim( x ); % 使用函数,根据前一步x的数值去计算对应的f和导数fprim,然后更新
  26.     iter = iter + 1;
  27.     fprintf('Iteration %d: x=%.20f, err=%.20f\n',iter, x, abs( f( x ) ) );
  28.     pause(.1);
  29. end
复制代码
回复此楼

wikys 发表于 2021-11-20 08:10:52
TouAkira 发表于 2021-11-20 07:58
仔细看循环,里面每轮迭代的更新,根本就没有使用到 f / f‘ 在当前x值时的取值,也就是说,在最关键的地方 ...

啊!!豁然开朗,感谢大佬的愿意帮我解答,不好意思之前没发现这么笨的错误!!!确实是看了几个视频就觉得自己好像懂了,开始套用模版刻舟求剑。。

那请问一下之前我循环的都是些什么东西呢?有点好奇。

还有就是我注意到您用了diff命令求导函数,这个和寻常的求导有什么不同吗?
我们学长说这个算式比较特殊,建议我们用手算...人工求导,他提到matlab得用特别的方式才行,所以我有点好奇~
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表