0
$\begingroup$

I have an iterative sum from $k=0$ to $k=n$ where the resulting sum is a polynomial of degree $n$.

I want to find the numerical root of this polynomial using FindRoot, starting from $x_0$ where the root I'm interested in is converging to some value. However I only know this $x_0$ for small $n \sim 60$.

Then FindRoot will give me a more precise $x_0$ and I can use this to find the root for $n=65$. I then use this root as my $x_0$ for $n=70$ etc. etc.

However I can only really increase $n$ by $5$ each time as I don't know where the root is converging and if I increase $n$ by more, FindRoot starts finding other roots which I know aren't the value I'm looking for.

How I can set it up so that I can run to $n=300$, for example, where FindRoot will use the root found $5$ iterations ago as its $x_0$?

Here's a sample:

max=90;
f[[0]]=1;

For[n=1,n<=max,n=n+1,
f[[n]]=Sum[f[[k]]Coefficient[S,y^k],{k,0,n}]];

MySum=Sum[f[[i]],{i,0,max}]

FindRoot[MySum ,{x,x_0},WorkingPrecision->7]

where $S$ is some polynomial $S(x,y)$

Improve this question
$\endgroup$
20
  • 1
    $\begingroup$ Please add working code to make this question easier to understand and answer and more useful for other visitors. $\endgroup$
    – Yves Klett
    Aug 20, 2014 at 9:46
  • $\begingroup$ @YvesKlett I would have. But my code are far too specific, tailored and long for it to be helpful for others. It's using about 20 predefined functions and impossible to understand without the whole notebook. Apologies $\endgroup$
    – Phibert
    Aug 20, 2014 at 9:57
  • 1
    $\begingroup$ Probably you can try then to define a minimum working example that includes the features you need. $\endgroup$
    – rhermans
    Aug 20, 2014 at 10:07
  • $\begingroup$ You should definitely go a for a MWE then. Without any code, this will likely be closed. Also, having to make up sample code to fit your verbal explanations is not very attractive. $\endgroup$
    – Yves Klett
    Aug 20, 2014 at 10:50
  • 1
    $\begingroup$ What is x_0? I don't think that's what you want. $\endgroup$
    – RunnyKine
    Aug 20, 2014 at 11:43

2 Answers 2

Reset to default
1
$\begingroup$

here is an example using a prior solution as the starting value for the next step:

polyn[x_?NumericQ, n_] := Normal@Series[ Sin[y], {y, 0, n}] /. y -> x
sol[7] = x /. FindRoot[ polyn[x, 7], {x, 3}]
Do[ sol[n] = x /. FindRoot[ polyn[x, n], {x, sol[n - 1]}], {n, 8, 20} ]

DiscretePlot[sol[n] - Pi, {n, 7, 20}, PlotRange -> All]

enter image description here

Improve this answer
$\endgroup$
1
  • $\begingroup$ This is exactly what I was looking for. $\endgroup$
    – Phibert
    Aug 21, 2014 at 10:45
1
$\begingroup$

This runs without errors. One definitely should not use the Head of an expression (e.g. f[[0]]) as a numeric variable in an array as the OP does, but I didn't want to rewrite everything. One normally uses f[[k+1]] to correspond to degree k and program the off-by-one indices throughout. Then the other Mathematica functions that work so efficiently on List can be used in your project.

max = 90;
f = ConstantArray[0, max];
S = x^2 + 3 x y + x y^3;
f[[0]] = 1;

For[n = 1, n <= max, n = n + 1, 
  f[[n]] = Sum[f[[k]] Coefficient[S, y, k], {k, 0, n}]];

MySum = Sum[f[[i]], {i, 0, max}]

FindRoot[MySum, {x, -1}, WorkingPrecision -> 7]

(*
  1 + 90 x^2 + 267 x^3 + 87 x (x^2 + 3 x^3)
  {x -> -1.011827}
*)

Finding a good starting point for FindRoot might be a problem. Since MySum is a polynomial, one can also use NSolve:

NSolve[MySum, x]
(*
  {{x -> -1.01183}, {x -> -0.380488},
   {x -> 0.0179959 - 0.0981234 I}, {x -> 0.0179959 + 0.0981234 I}}
*)
Improve this answer
$\endgroup$
5
  • $\begingroup$ Thank you very much for working with my slightly void question. So yes I use NSolve to find a good starting point (I should note: I don't use Nsolve for very large max because it breaks down as my polynomial is too large. Hence FindRoot seems to be my only option.) However I can't see how this code updates $x_0$ as $n$ runs from $0$ to $max$. If I use the "starting point" the whole way, the root quickly gets too far from this value and FindRoot starts finding other roots which I'm not interested in. $\endgroup$
    – Phibert
    Aug 20, 2014 at 14:51
  • $\begingroup$ @user13223423 You provided no code for $x_0$. I have no idea what you want to do with it. The way you use it, it is just a starting point for FindRoot. You'll have to figure out how to get a close enough starting point to the root you seek. By the way, you should write it in Mathematica as x0. The expression x_0 is a Pattern, not a number. $\endgroup$
    – Michael E2
    Aug 20, 2014 at 14:59
  • $\begingroup$ For each max, my polynomial has a different root. As I increase max (manually), the root starts moving until soon the root it is too far away from my starting point for FindRoot to locate it. My question is: how can I can set this up so that the starting point for FindRoot changes each time $n$ increases from $0$ to max, and uses the root found from the same calculation but max-1. This way, my "starting point" will always be close enough to the root I want for FindRoot to find it. But sorry, this seems to have wasted your time $\endgroup$
    – Phibert
    Aug 20, 2014 at 15:11
  • $\begingroup$ @user13223423 If you can say what to change x0 to each time, I can write code to do that. I don't have time to solve the problem of how to determine the starting point, and one might need to know all the details of the problem to do that. You seem to know something about how it moves. Can you not use that to get somewhat close? Or can you always find two numbers, x1, x2 that bracket the root (for which the polynomial is pos. and neg.). If so, look at the FindRoot help. E.g. FindRoot[f, {x, x0, x1, x2}]. $\endgroup$
    – Michael E2
    Aug 20, 2014 at 15:23
  • $\begingroup$ I'll give that a go, thanks for the help $\endgroup$
    – Phibert
    Aug 20, 2014 at 16:29

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.