# Having trouble translating Matlab fsolve into Mathematica

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? Jun 30, 2015 at 16:07
• @J. M. The description of the original problem and code is here, jingplusplus.com/2015/06/12/mathpuzzleslotteryproblem
– Nick
Jun 30, 2015 at 16:10

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} *)


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.

• 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, 2015 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. Jun 30, 2015 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 1, 2015 at 3: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. Jul 1, 2015 at 3: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 1, 2015 at 4:47

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}}} *)


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.