如何在 Octave 中使用 for-end 和 if-else 计算值？

问题描述

我需要在以下等式中计算变量 Nfj：

ea_M = (sf/E)*(2*Nfj)^b + ef*(2*Nfj)^c

所以它等于 ea 的 4 个不同值：5.0900e-04、4.3626e-04、3.6358e-04 和 2.9084e-04。结果应该是 Nfj 的 4 个值，应该存储在 results 中。我四舍五入是因为我认为不可能计算出完全相等的 ea 值，如果不四舍五入，则相等总是错误的。

我为它编写了这段代码，但它没有按预期工作：脚本运行了很长时间没有任何结果。如何修复它才能正常工作？

sf = 882.07;
ef = 0.59;
b = -0.102969;
c = -0.58;
E = 210000;
ea = [5.0900e-04; 4.3626e-04; 3.6358e-04; 2.9084e-04]

for pos = 1:length(ea)
  for Nfj = 1e3:10:1e12
    ea_M = (sf/E)*(2*Nfj)^b + ef*(2*Nfj)^c;
    if  round(ea_M * 10^5)/10^5 == round(ea(pos) * 10^5)/10^5;
      disp(ea_M)
      disp(Nfj)
      results(pos) = Nfj;
    end
  end
end

解决方法

您正在尝试每 10^th 个值，从 10³ 到 10¹²：有 10¹¹值这里试试！这当然很多。您将永远寻找，并且您可能会跳过使等式成立的实际值。

如果您无法手动求解方程，可以使用数值求解器。我们先打印函数：

sf = 882.07;
ef = 0.59;
b = -0.102969;
c = -0.58;
E = 210000;

ea_M = @(Nfj) (sf/E)*(2*Nfj).^b + ef*(2*Nfj).^c;

Nfj = logspace(3,12,1000);
plot(Nfj,ea_M(Nfj))
set(gca,'xscale','log')

看起来这个函数非常单调，您要查找的四个值位于您正在搜索的 10³ 到 10¹² 的区间内。要找到它等于您的值之一的位置，我们可以减去该值并找到它等于零的位置。如果您从函数大于零的点和函数小于零的点开始，您可以非常快速地缩小搜索范围。您每次将间隔减半，保留包含过零的一半间隔。 fzero 函数就是这样做的。

ea = [5.0900e-04; 4.3626e-04; 3.6358e-04; 2.9084e-04];
results = zeros(size(ea));
for pos = 1:numel(ea)
   results(pos) = fzero(@(Nfj) ea_M(Nfj) - ea(pos),[1e3,1e12]);
end
results

在 MATLAB 中，此代码在几分之一秒内运行并输出：

results =

   1.0e+10 *

    0.0429
    0.1849
    1.0627
    9.1919

并且 ea_M(results) - ea 近似为零。

关于我在这里发布的代码的一些说明：

ea_M 被定义为匿名函数。这样可以更轻松地重用表达式，而不是一遍又一遍地编写它。我将 ^ 替换为 .^ 以允许此函数一次对 Nfj 值的数组进行计算，而不仅仅是单个值。这是 fzero 调用所必需的。
我在对数刻度上绘制了该函数，因为该函数需要这样做。并非所有函数都如此。
我预先分配了 results 数组，您应该避免在循环中增加数组的大小。

octave octave