Am I using a wrong numerical method?

Question

This is the code:

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')

If I remove 'AllRealRoots' option it works fast and finds a solution, but when I enable the option Matlab does not finish for an hour. Am I using a wrong numerical method?

Answer 1

First, straight from the documentation for numeric::solve:

If eqs is a non-polynomial/non-rational equation or a set or list containing such an equation, then the equations and the appropriate optional arguments are passed to the numerical solver numeric::fsolve.

So, as your equation f is non-polynomial, you should probably call numeric::fsolve directly. However, even with the 'MultiSolutions' it fails to return more than one root over your range (A bug perhaps? – I'm using R2013b). A workaround is to call numeric::realroots to get bounds on each of the district real roots in your range and then solve for them separately:

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

The resultant solution vector, T, is double-precision (allocate it as sym or vpa you want higher precision):

T =

  -0.663159371123072
   0.034848320470578
   0.999047064621451
   2.000000000000000
   2.695929753727520
   3.933983894260340
   4.405822476913172
   5.868112290810963
   6.108685019679461

You may be able to remove the for loop if you can figure out how to cleanly pass in the output of 'numeric::realroots' to 'numeric::fsolve' in one go (it's doable, but may require converting stuf to strings unless you're clever).

Another (possibly even faster) approach is to switch to using the numeric (floating-point) function fzero for the second half after you bound all of the roots:

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

I checked and, for your problem here and using the default tolerances, the resultant T is within a few times machine epsilon (eps) of the numeric::fsolve' solution.