我使用了错误的数值方法吗？答案

【问题标题】：Am I using a wrong numerical method?我使用了错误的数值方法吗？
【发布时间】：2014-02-14 17:45:41
【问题描述】：

这是代码：

f = dsolve('D3y+12*Dy+y = 0 ,y(2) = 1 ,Dy(2) = 1, D2y(2) = -1');
feval(symengine, 'numeric::solve',strcat(char(f),'=1'),'t=-4..16','AllRealRoots')

如果我删除 'AllRealRoots' 选项，它可以快速运行并找到解决方案，但是当我启用该选项时，Matlab 不会在一个小时内完成。我是否使用了错误的数值方法？

【问题讨论】：

标签： performance matlab numeric symbolic-math

【解决方案1】：

首先，直接来自numeric::solve 的文档：

如果 eqs 是非多项式/非有理方程或包含此类方程的集合或列表，则将方程和适当的可选参数传递给数值求解器 numeric::fsolve。

因此，由于您的方程 f 是非多项式的，您可能应该直接调用 numeric::fsolve。但是，即使使用'MultiSolutions'，它也无法在您的范围内返回一个以上的根（也许是一个错误？-我正在使用 R2013b）。一种解决方法是调用numeric::realroots 以获取您范围内每个区域实际根的界限，然后分别解决它们：

f = dsolve('D3y+12*Dy+y = 0 ,y(2) = 1 ,Dy(2) = 1, D2y(2) = -1');
r = feval(symengine, 'numeric::realroots', f==1, 't = -4 .. 16');

num_roots = numel(r);
T = zeros(num_roots,1); % Wrap in sym or vpa for higher precision output
syms t;
for i = 1:num_roots
    bnds = r(i);
    ri = feval(symengine, '_range', bnds(1), bnds(2));
    s = feval(symengine, 'numeric::fsolve', f==1, t==ri);
    T(i) = feval(symengine, 'rhs', s(1));
end

得到的解向量 T 是双精度的（分配为 sym 或 vpa 你想要更高的精度）：

T =

  -0.663159371123072
   0.034848320470578
   0.999047064621451
   2.000000000000000
   2.695929753727520
   3.933983894260340
   4.405822476913172
   5.868112290810963
   6.108685019679461

如果您能弄清楚如何一次性将'numeric::realroots' 的输出干净地传递给'numeric::fsolve'，您也许可以删除for 循环（这是可行的，但可能需要将stuf 转换为字符串，除非您'很聪明）。

另一种（可能更快）方法是在绑定所有根后切换到使用数字（浮点）函数fzero 处理下半部分：

f = dsolve('D3y+12*Dy+y = 0 ,y(2) = 1 ,Dy(2) = 1, D2y(2) = -1');
r = feval(symengine, 'numeric::realroots', f==1, 't = -4 .. 16');

num_roots = numel(r);
T = zeros(num_roots,1);
g = matlabFunction(f-1); % Create anonymous function from f
for i = 1:num_roots
    bnds = double(r(i));
    T(i) = fzero(g,bnds);
end

我检查了，对于您在此处的问题并使用默认容差，结果 T 是 numeric::fsolve' 解决方案的几倍机器 epsilon (eps)。

【讨论】：