Message Boards Message Boards

0
|
5983 Views
|
10 Replies
|
3 Total Likes
View groups...
Share
Share this post:

How to add step monitor to Nsolve?

Posted 12 years ago
This is my piece of code
For[x = 1/2; f = k*x*(1 - x); i = 1, i < (2^4), i++, g = k*f*(1 - f);
f = g; Print[i][i];]
rts = NSolve[(g - 1/2) == 0, {k}, WorkingPrecision -> 10]
I wish to modify this code to
  1. To continuously output rts matrix as it solves the equation
  2. Show a progress bar to how how long it may take
  3. If possible to show roots around 3.5  >> Note : because of the high degree of  the recursive equation FindRoot is not working
POSTED BY: ashwin swarup
10 Replies
Posted 12 years ago
Hmm.. i knew it would be larger ..but not that large .

I have an estimate of where the root i require should be.

Plotting is not really that important
I need it to find the next root closest to 3.5

But i guess its expecting too much from the program .
Ill see what i can do

Thank you for your help learnt something new emoticon
POSTED BY: ashwin swarup
Posted 12 years ago
Each iteration of your For loop approximately doubles the number of terms in your f expression.

LeafCount will count the number of nodes to make up an expression. Leafcount[3+4*x]==5.
When you had i<2*4 LeafCount was 262142. Now with i<2^5 LeafCount is 17179869182.
You can verify this by changing your Print to Print[i, " ", LeafCount]].

Plot will need to evaluate f many times to create the plot. Consider evaluating 17179869182 items many times.
I would estimate plotting i<2^5 will take 2^16 times longer than plotting i<2^4. This may take 1000 hours to Plot!
But if Mathematica uses all available memory and must start storing results on your hard drive swap space
then the calculation can take 100 or 1000 times longer to complete.
I suspect Mathematica or the computer may crash long before it finishes the Plot.
You can estimate how long it will take by timing i<12, then i<13, then i<14, then i<15 and drawing a graph of the times.

You might save a small amount of time and by changing
g = k*f*(1 - f); f = g;
...
Plot[g....
...
FindRoot[g...

to

f = k*f*(1 - f);

...
Plot[f...
...
FindRoot[f...

That may save memory by not needing to have both f and g in memory.
But this will probably not help enough.

You can use Put to store the value of g in a file. Later you can use Get to read the value of g from the file.
The help system shows you how to use these.

Sometimes results that are interesting are worth waiting for, if you can estimate they will be done in a week or a month.

Using the fastest single core computer with as much memory as it can use may make i<2^5 go many many times faster.
Memory is inexpensive. I suggest 32 gigabytes for working on large or hard problems in Mathematica.
POSTED BY: Bill Simpson
Posted 12 years ago
Another thing is once g is evaluated is it possible to permanently store the function as a file so that it doesnt have to keep calculating g every time i have to run the program ?

Also i converted to numeric function as suggested by Youngjoo Chung  . This too didnt help the cause
POSTED BY: ashwin swarup
Posted 12 years ago
My machine runs 8.1 mathematica
It has 4 GB RAM
 
 x = 1/2;
 f = k*x*(1 - x);
 For[i = 1, i < 2^5, i++, g = k*f*(1 - f); f = g; Print[i];];
 Print["Finished calculating g, starting plot"];
 Plot[g - 1/2, {k, 3.560, 3.575}, WorkingPrecision -> 50,
  PlotRange -> All]
 Print["Finished plot, starting findroot"];
 FindRoot[g - 1/2 == 0, {k, 3567/1000, 3566/1000, 3569/1000},
WorkingPrecision -> 50, StepMonitor :> Print[k]][/i]
As you notice the equation is 2^5 >
Also te thing is getting stuck at after calculating g i.e. at the plot.
I reduced wroking precision to 50 and then 20 with no result .
...Wonder what the way out could be.
POSTED BY: Updating Name
Hi,

FindRoot takes a long time because it is trying to do symbolic evaluation with g, a huge polynomial. You can speed up calculation by making a function definition with restriction. Try the following:
gfun[k_?NumericQ] := g;
Plot[gfun[k] - 1/2, {k, 0, 4}, AxesLabel -> {k, "g"}]
{FindRoot[(gfun[k] - 1/2) == 0, {k, 1.9, 2.1}],
FindRoot[(gfun[k] - 1/2) == 0, {k, 3.2, 3.25}],
FindRoot[(gfun[k] - 1/2) == 0, {k, 3.45, 3.5}]}

I got the numbers for k in FindRoot from the plot. g is extremely oscillatory about zero beyond k=3.5. I got the result below.
In[563]:= {FindRoot[(gfun[k] - 1/2) == 0, {k, 1.9, 2.1}],
  FindRoot[(gfun[k] - 1/2) == 0, {k, 3.2, 3.25}],
  FindRoot[(gfun[k] - 1/2) == 0, {k, 3.45, 3.5}]} // Timing

Out[563]= {0.937500, {{k -> 2.}, {k -> 3.23607}, {k -> 3.49856}}}

Youngjoo Chung
POSTED BY: Youngjoo Chung
Posted 12 years ago
Is there a way of knowing where the above code is getting stuck

>> If it is getting stuck in the first place
Or some sort of indication that i have to wait it out and not restart mathematica
>> I thought step monitor did that ..but doesnt seem to be doing so

tyvm
POSTED BY: ashwin swarup
Posted 12 years ago
I modified the code
For[x = 1/2; f = k*x*(1 - x); i = 1, i < (2^5), i++, g = k*f*(1 - f);
f = g; Print[i][i];]
Plot[g - 1/2, {k, 3.3, 3.7}, WorkingPrecision -> 100, PlotRange -> All]
FindRoot[g - 1/2 == 0, {k, 350/100, 348/100, 352/100},
WorkingPrecision -> 100, StepMonitor :> Print[k]][/i]

Its getting stuck and not able to plot the  function
Is the function too big for Mathematica to handle ?
POSTED BY: ashwin swarup
Posted 12 years ago
Hi can u explain to me what "350/100, 348/100, 352/100" are for ? what do they do ?

From http://reference.wolfram.com/mathematica/ref/FindRoot.html
click on "Details and Options"  to see hidden under that it says
FindRoot[lhs==rhs {x,xstart, xmin, xmax}]
 searches for a solution, stopping the search if x ever gets outside the range xmin to xmax.

Suggestion: Perhaps that line should be moved above Details and Options to be more visible because it is very often not found by users and is commonly needed.

When I look at the plot I see your function is smooth between 3.48 and 3.52 with a root between those.
This should only be used when you can see your function is well behaved between xmin and xmax and a solution will be simply found. 
I use this to force FindRoot to only search very near to 3.5 and not accidentally begin looking for roots beyond 3.6 where it would become hopelessly confused and take forever.

If you replace 348/100 and 352/100 by 3.48 and 3.52 you will see Findroot begins with much less than 100 digits of precision.
348/100 and 352/100 have infinite precision. You could also use 3.48`100 and 3.52`100 to give 100 digits of precision. 

Its getting stuck and not able to plot the function
Is the function too big for Mathematica to handle ?


Perhaps it is busy calculating g. Perhaps it is busy calculating the plot.
Try this and see if it helps you determine where the time is being spent.
x = 1/2;
f = k*x*(1 - x);
For[i = 1, i < 2^4, i++, g = k*f*(1 - f); f = g];
Print["Finished calculating g, starting plot"];
Plot[g - 1/2, {k, 3.3, 3.7}, WorkingPrecision -> 100, PlotRange -> All]
Print["Finished plot, starting findroot"];
FindRoot[g - 1/2 == 0, {k, 350/100, 348/100, 352/100}, WorkingPrecision -> 100, StepMonitor :> Print[k]]

How fast is your computer? How much memory does your computer have? What version of Mathematica do you have?

If the time is taken by the plot then try reducing WorkingPrecision to 50 or even 20.
Very near 3.5, less than 3.6 this will not change the Plot substantially.
Beyond 3.6 more precision is required to see more accurate results.
If the Plot is still very slow try Plot from 0 to 3.55 and see if the plot above 3.6 is taking the time.

If the time is taken in FindRoot then again try reducing WorkingPrecision to 50 or even 20.
Then increase the WorkingPrecision gradually and see if the answer stays the same or if it changes.
If the answer changes then that is a warning for you to spend more time working on an accurate solution.
POSTED BY: Bill Simpson
Posted 12 years ago
Hi can u explain to me what  "350/100, 348/100, 352/100" are for ? what do they do ?
POSTED BY: ashwin swarup
Posted 12 years ago
x=1/2;
f=k*x*(1-x);
For[i=1, i<2^4, i++,
  g=k*f*(1-f);
  f=g
  ];
Plot[g-1/2, {k, 3.3, 3.7}, WorkingPrecision->100, PlotRange->All]
FindRoot[g-1/2==0, {k, 350/100, 348/100, 352/100}, WorkingPrecision->100, StepMonitor :> Print[k]]
POSTED BY: Bill Simpson
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract