matlab中的严重错误解决 - 我该怎么办？

Question

这篇文章在Solve中更新我的问题找到错误的解决方案？

鉴于x中的"简单"功能:

f = (x + 3/10)^(1/2) - (9*(x + 3/10)^5)/5 - (x + 3/10)^6 + (x - 1)/(2*(x + 3/10)^(1/2));

通过电话找到零

solve(f,x)

这产生3个零:

 ans =

 0.42846617518653978966562924618638
 0.15249587894102346284238111155954
 0.12068186494007759990714181154349

简单地看一下图表显示第三个根是无意义的:

我有一个严重的问题,因为我需要从上面的向量中检索最小的零.调用min(ans)返回错误的零.我该怎么做才能解决？

Answer 1

这是一个非多项式方程,它可能会回退到数值解算器(非符号).所以可能存在数字错误,或者数字算法可能会卡住并报告错误的解决方案,我不确定......

您可以做的是将解决方案替换回等式,并拒绝高于某个指定阈值的解决方案:

% define function
syms x real
syms f(x)
xx = x+3/10;
f(x) = sqrt(xx) - 9/5*xx^5 - xx^6 + (x - 1)/(2*sqrt(xx));
pretty(f)

% find roots
sol = solve(f==0, x, 'Real',true)

% filter bad solutions
err = subs(f, x, sol)
sol = sol(abs(err) < 1e-9);  % this test removes the 2nd solution

% plot
h = ezplot(f, [0.1 0.5]);
line(xlim(), [0 0], 'Color','r', 'LineStyle',':')
xlabel('x'), ylabel('f(x)')

% programmatically insert data tooltips
xd = get(h, 'XData'); yd = get(h ,'YData');
[~,idx] = min(abs(bsxfun(@minus, xd, double(sol))), [], 2);
dcm = datacursormode(gcf);
pos = [xd(idx) ; yd(idx)].';
for i=1:numel(idx)
    dtip = createDatatip(dcm, h);
    set(get(dtip,'DataCursor'), 'DataIndex',idx(i), 'TargetPoint',pos(i,:))
    set(dtip, 'Position',pos(i,:))
end

我们只剩下两个期望的解决方案(一个被我们的测试拒绝):

                   /      3 \5
                 9 | x + -- |
    /      3 \     \     10 /    /      3 \6         x - 1
sqrt| x + -- | - ------------- - | x + -- |  + ----------------
    \     10 /         5         \     10 /          /      3 \
                                               2 sqrt| x + -- |
                                                     \     10 /
sol =
 0.42846617518653978966562924618638
 0.12068186494007759990714181154349       % <== this one is dropped
 0.15249587894102346284238111155954

err(x) =
 -9.1835496157991211560057541970488e-41
   -0.058517436737550288309001512815475   % <==
  1.8367099231598242312011508394098e-40

我也尝试过使用MATLAB的数值求解器,它们能够找到合理起点的两个解决方案:

• fzero函数(无导数的根寻找算法)
• fsolve功能(来自优化工具箱)

(见相关问题)

fcn = matlabFunction(f);  % convert symbolic f to a regular function handle
opts = optimset('Display','off', 'TolFun',1e-9, 'TolX',1e-6);

% FZERO
sol2(1) = fzero(fcn, 0.1, opts);
sol2(2) = fzero(fcn, 0.5, opts);
disp(sol2)    % [0.1525, 0.4285]

% FSOLVE
sol3(1) = fsolve(fcn, 0.0, opts);
sol3(2) = fsolve(fcn, 1.0, opts);
disp(sol3)    % [0.1525, 0.4285]

为了比较,我尝试直接将方程式解释为MuPAD以及Mathematica.

Mathematica 9.0

MuPAD(R2014a)

当然我们总是可以直接从MATLAB中调用MuPAD:

% f is the same symbolic function we've built above

>> sol = evalin(symengine, ['numeric::solve(' char(f) ' = 0, x, AllRealRoots)'])
sol =
[ 0.15249587894102346284238111155954, 0.42846617518653978966562924618638]

上述调用相当于搜索整个范围x = -infinity .. infinity(可能很慢!).我们应该numeric::solve尽可能提供特定的搜索范围:

>> sol = feval(symengine, 'numeric::solve', f==0, 'x = 0 .. 1', 'AllRealRoots')
sol =
[ 0.15249587894102346284238111155954, 0.42846617518653978966562924618638]