1
$\begingroup$

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)]
$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – bbgodfrey Jun 30 '15 at 16:02
  • $\begingroup$ @bbgodfrey myfun14is included in the code parts $\endgroup$ – Nick Jun 30 '15 at 16:04
  • $\begingroup$ What exactly are you trying to do to the companion matrix of a polynomial? $\endgroup$ – J. M. is away Jun 30 '15 at 16:07
  • $\begingroup$ @Guesswhoitis. The description of the original problem and code is here, jingplusplus.com/2015/06/12/mathpuzzleslotteryproblem $\endgroup$ – Nick Jun 30 '15 at 16:10
4
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Nick Jun 30 '15 at 17:02
  • $\begingroup$ @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. $\endgroup$ – bbgodfrey Jun 30 '15 at 17:57
  • $\begingroup$ 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. $\endgroup$ – Nick Jul 1 '15 at 3:16
  • $\begingroup$ @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. $\endgroup$ – bbgodfrey Jul 1 '15 at 3:55
  • $\begingroup$ 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... $\endgroup$ – Nick Jul 1 '15 at 4:47
0
$\begingroup$

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.