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

=================

  

 

Welcome to Mathematica.SE! I hope you will become a regular contributor. To get started, 1) take the introductory Tour now, 2) when you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge, 3) remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign, and 4) give help too, by answering questions in your areas of expertise.
– bbgodfrey
Jun 30 ’15 at 16:02

  

 

@bbgodfrey myfun14is included in the code parts
– Nick
Jun 30 ’15 at 16:04

  

 

What exactly are you trying to do to the companion matrix of a polynomial?
– J. M.♦
Jun 30 ’15 at 16:07

  

 

@Guesswhoitis. The description of the original problem and code is here, jingplusplus.com/2015/06/12/mathpuzzleslotteryproblem
– Nick
Jun 30 ’15 at 16:10

=================

2 Answers
2

=================

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.

  

 

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

  

 

Thanks! I will do some check.
– Nick
Jul 1 ’15 at 3:06

  

 

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 ’15 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.
– bbgodfrey
Jul 1 ’15 at 3:55

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.