Having trouble translating Matlab fsolve into Mathematica

Question

I searched the site and found that someone said FindRoot is equivalent to fsolve, however, without proper options, FindRoot is going to run forever. In my case, I am trying to do the numerical approximations, but it seems that Mathematica is doing a lot of symbolic calculations.

Here is the original Matlab code I want to translate:

x0 = 0.06*ones(14,1);
options = optimoptions('fsolve','Display','iter','TolFun',1e-60,'TolX',1e-60);
[x1,fval14] = fsolve(@myfun14,x0,options);

And my attempt is:

x1=FindRoot[myfun14[Array[x, 14]] == 0, Transpose[{Array[x, 14], ConstantArray[0.06, 14]}]]

Thanks!

The full original matlab code and the mathematica code I translated are as follows for reference:

Matlab code:

function F14 = myfun14(x)
  b = [1/2;1/4;0;1/8;0;0;0;1/16;0;0;0;0;0;0];
  A = [0,0,0,0,0,0,0,0,0,0,0,0,0,x(1);
  1,0,0,0,0,0,0,0,0,0,0,0,0,x(2);
  0,1,0,0,0,0,0,0,0,0,0,0,0,x(3);
  0,0,1,0,0,0,0,0,0,0,0,0,0,x(4);
  0,0,0,1,0,0,0,0,0,0,0,0,0,x(5);
  0,0,0,0,1,0,0,0,0,0,0,0,0,x(6);
  0,0,0,0,0,1,0,0,0,0,0,0,0,x(7);
  0,0,0,0,0,0,1,0,0,0,0,0,0,x(8);
  0,0,0,0,0,0,0,1,0,0,0,0,0,x(9);
  0,0,0,0,0,0,0,0,1,0,0,0,0,x(10);
  0,0,0,0,0,0,0,0,0,1,0,0,0,x(11);
  0,0,0,0,0,0,0,0,0,0,1,0,0,x(12);
  0,0,0,0,0,0,0,0,0,0,0,1,0,x(13);
  0,0,0,0,0,0,0,0,0,0,0,0,1,x(14)];
      B = A^(16-15);
      for n=5:50
      B = B + A^(2^n-15)/2^(n-4);
      end
      F14 = (x - b - 1/32 * B * x);


x0 = 0.06*ones(14,1);
options = optimoptions('fsolve','Display','iter','TolFun',1e-60,'TolX',1e-60);
[x1,fval14] = fsolve(@myfun14,x0,options);
[x2,fval14] = fsolve(@myfun14,x1,options);
[x3,fval14] = fsolve(@myfun14,x2,options);
[x4,fval14] = fsolve(@myfun14,x3,options);
[x5,fval14] = fsolve(@myfun14,x4,options);
[x,fval14] = fsolve(@myfun14,x5,options);
for t=1:200
[x,fval14] = fsolve(@myfun14,x,options);
end

A = [0,0,0,0,0,0,0,0,0,0,0,0,0,x(1);
1,0,0,0,0,0,0,0,0,0,0,0,0,x(2);
0,1,0,0,0,0,0,0,0,0,0,0,0,x(3);
0,0,1,0,0,0,0,0,0,0,0,0,0,x(4);
0,0,0,1,0,0,0,0,0,0,0,0,0,x(5);
0,0,0,0,1,0,0,0,0,0,0,0,0,x(6);
0,0,0,0,0,1,0,0,0,0,0,0,0,x(7);
0,0,0,0,0,0,1,0,0,0,0,0,0,x(8);
0,0,0,0,0,0,0,1,0,0,0,0,0,x(9);
0,0,0,0,0,0,0,0,1,0,0,0,0,x(10);
0,0,0,0,0,0,0,0,0,1,0,0,0,x(11);
0,0,0,0,0,0,0,0,0,0,1,0,0,x(12);
0,0,0,0,0,0,0,0,0,0,0,1,0,x(13);
0,0,0,0,0,0,0,0,0,0,0,0,1,x(14)];
xfinal = (A)^((10^9)-15)*x;

Mathematica code for the function definition part:

myfun14[list_] := Module[{b, A, B},
b = Transpose@{{1/2., 1/4., 0., 1/8., 0, 0, 0, 1/16., 0, 0, 0, 0, 0,
   0}};
A = {{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, list[[1]]},
{1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, list[[2]]},
{0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, list[[3]]},
{0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, list[[4]]},
{0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, list[[5]]},
{0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, list[[6]]},
{0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, list[[7]]},
{0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, list[[8]]},
{0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, list[[9]]},
{0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, list[[10]]},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, list[[11]]},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, list[[12]]},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, list[[13]]},
{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, list[[14]]}};
B = A + Sum[MatrixPower[A, (2^n - 15)]/2.^(n - 4), {n, 5, 50}];
(list - b - 1/32.*B*list)]

What exactly are you trying to do to the companion matrix of a polynomial? — J. M.'s missing motivation, Jun 30 '15 at 16:07
@J. M. The description of the original problem and code is here, http://jingplusplus.com/2015/06/12/mathpuzzleslotteryproblem/ — Nick, Jun 30 '15 at 16:10

bbgodfrey · Accepted Answer · 2015-07-01T14:31:47.333

7

Changing B*list to B.list seems to solve the problem. Also, I recommend dropping unnecessary decimal points from myfun14. Doing so gives

(* {x[1] -> 0.5, x[2] -> 0.273438, x[3] -> 0.0128174, x[4] -> 0.125601, 
    x[5] -> 0.00588754, x[6] -> 0.000275978, x[7] -> 0.0000129365, 
    x[8] -> 0.0625006, x[9] -> 0.00292972, x[10] -> 0.00013733, 
    x[11] -> 6.43736*10^-6, x[12] -> 3.01751*10^-7, 
    x[13] -> 1.41446*10^-8, x[14] -> 6.63028*10^-10} *)

Addendum

The recently updated code in the Question truly takes forever to run, because FindRoot first attempts to evaluate symbolically the gigantic expression created by myfun14[Array[x, 14]]. This can be circumvented by using the option, Evaluated -> False.

FindRoot[myfun14[Array[x, 14]], Transpose[{Array[x, 14], ConstantArray[6/100, 14]}], 
 Evaluated -> False]
(* {x[1] -> 0.500634, x[2] -> 0.266567, x[3] -> 0.00914348, 
    x[4] -> 0.130162, x[5] -> 0.00916464, x[6] -> 0.00235785, 
    x[7] -> 0.00297105, x[8] -> 0.065152, x[9] -> 0.00401983, 
    x[10] -> 0.00206912, x[11] -> 0.00273191, x[12] -> 0.00231037, 
    x[13] -> 0.00135524, x[14] -> 0.00136153} *)

Runtime on my PC is just less than one second.

edited Jul 01 '15 at 14:31

answered Jun 30 '15 at 16:28

bbgodfrey

61,439
17
89
156

Thank you, your method worked. I still have two things not clear : 1. I added those decimal points to let mathematica do numerical computation instead of keeping the rational numbers, but I have made things worse; 2. Although the program seems to be identical with the matlab one, this step gives different and wrong answers compared to matlab, where the correct answer should be approximately {0.5006, 0.2666, 0.0091, 0.1302, 0.0092, 0.0024, 0.0030, 0.0652, 0.0040, 0.0021, 0.0027, 0.0023, 0.0014, 0.0014}. Can you explain a little bit? – Nick Jun 30 '15 at 17:02
@Nick 1) The problem appears to have been MatrixPower[A, (16 - 15.)]. MatrixPower expects an integer for its second argument, and 15. is a real number. 2} I tried different initial guesses and a much higher WorkingPrecision for FindRoot, but the answer did not change. My guess is that there is a typo somewhere. – bbgodfrey Jun 30 '15 at 17:57
I have located the typo, I forgot to change the ^ into MatrixPower in the definition of myFun14. But after changing this, the question returned to its original state: the function again becomes evaluating forever:( I have updated the mathematica code. – Nick Jul 01 '15 at 03:16
@Nick Just a guess: replace 2.^(n - 4) by 2^(n - 4). In general, use rationale numbers, unless there is a good reason not to. – bbgodfrey Jul 01 '15 at 03:55
Unfortunately it does not work this time. My guess is that MatrixPower is trying to evaluate the power symbolically and then guess numerically using FindRoot, while matlab is doing both numerically... – Nick Jul 01 '15 at 04:47
@Nick With an easy but less than obvious change, I get {x[1] -> 0.500634, x[2] -> 0.266567, x[3] -> 0.00914348, x[4] -> 0.130162, x[5] -> 0.00916464, x[6] -> 0.00235785, x[7] -> 0.00297105, x[8] -> 0.065152, x[9] -> 0.00401983, x[10] -> 0.00206912, x[11] -> 0.00273191, x[12] -> 0.00231037, x[13] -> 0.00135524, x[14] -> 0.00136153}} in 0.9 seconds. See edit to Answer. – bbgodfrey Jul 01 '15 at 14:26

score 2 · Answer 2 · answered Mar 22 '20 at 06:23

We can also make use of the undocumented syntax of FindRoot:

Clear@myfun14; 
myfun14[list_] := 
 Module[{b, A, B}, 
  b = Transpose@{{1/2, 1/4, 0, 1/8, 0, 0, 0, 1/16, 0, 0, 0, 0, 0, 0}};
  A = DiagonalMatrix[ConstantArray[1., 13], -1]; A[[All, -1]] = list;
  B = A + Sum[MatrixPower[A, (2^n - 15)]/2.^(n - 4), {n, 5, 50}];
  (list - b - 1/32 B.list)]

FindRoot[myfun14, {ConstantArray[6/100, 14]}] // AbsoluteTiming

(* {1.13405, {{0.500634, 0.266567, 0.00914348, 0.130162, 
               0.00916464, 0.00235785, 0.00297105, 0.065152, 0.00401983, 
               0.00206912, 0.00273191, 0.00231037, 0.00135524, 0.00136153}}} *)

score 1 · Answer 3 · answered Jul 04 '15 at 15:31

1

An alternative solution would be changing the myfun14[list_] to myfun14[list_ /; VectorQ[list, NumberQ]], the result can be obtained after a few seconds. Setting a high AccuracyGoal is also necessary to get the correct answer for the problem.

answered Jul 04 '15 at 15:31

vapor

7,911
2
22
55

Having trouble translating Matlab fsolve into Mathematica

3 Answers3