FindRoot inside Table

(MWE at the end of the post)

I need to solve a non-linear equation f(y;x1,x2,..,x5)f(y;x1,x2,..,x5)f(y;x_1,x_2,..,x_5) in one variable yyy and then compute 4 new output expressions, for over 60 different initial parameter inputs of xixix_i. The 4 output variables which are as follows:

gm(y,x1,..,x5)∀mgm(y,x1,..,x5)∀mg_m(y,x_1,..,x_5) \; \forall m

The main equation will solve for the variable y(i)y(i)y(i):

f(y(i),x1(i),x2(i),…,x5(i))=0∀i∈{1,60}f(y(i),x1(i),x2(i),…,x5(i))=0∀i∈{1,60}f(y(i),x_1(i),x_2(i),…,x_5(i))=0 \; \forall \; i \in \{1,60\}

Now my

f(y(i),x1(i),x2(i),…,x5(i))=A1(y,x1,..,x5)∗yk1+B1(y,x1,..,x5)∗yk2−c1(x1,..,x5)f(y(i),x1(i),x2(i),…,x5(i))=A1(y,x1,..,x5)∗yk1+B1(y,x1,..,x5)∗yk2−c1(x1,..,x5)f(y(i),x_1(i),x_2(i),…,x_5(i))=A1(y,x_1,..,x_5)*y^{k_1} + B1(y,x_1,..,x_5)*y^{k_2}-c1(x_1,..,x_5)

I define variables

A1=A1(y,x1,..,x5);
B1=B1(y,x1,..,x5);
c1=c1(x1,..,x5);
g1=g1(y,x1,..,x5);
g2=g2(y,x1,..,x5);
g3=g3(y,x1,..,x5);
g4=g4(x1,..,x5);

I have data for x1(i),..,x5(i)x1(i),..,x5(i)x_1(i),..,x_5(i) in a CSV file that I Import and Table straight into the variable names making them into lists.

datatemp = Import[“C:\\Documents\\2012U26G0.csv”];
j = Dimensions[datatemp][[1]]
kk=2
x1 = Table[datatemp[[i, 3]], {i, kk, j}]
x2 = Table[datatemp[[i, 4]], {i, kk, j}]
x3 = Table[datatemp[[i, 5]], {i, kk, j}]
x4 = Table[datatemp[[i, 6]], {i, kk, j}]
x5 = Table[datatemp[[i, 7]], {i, kk, j}]

I think this automatically makes my earlier defined formulae for A1A1A1,B1B1B1 and gmgmg_m into a list of formulae with the only unknown being yyy and it makes c1c1c1 into a list of constants since c1c1c1 was only dependent on xixix_i.

Now, what I would like to be able to do is the following.
Give some initial search point for my FindRoot for i=1i=1i=1

sol={1,1,1,1,120}

As you will see in a second, I only care about

sol[[5]]

Due to continuity, the roots move monotonically with iii so once I find one root, I can get a sense of where to look for the next one so I substitute the previous solution into the search for the next one. Also, in one shot I compute the 4 output variables I need. So when the Table is run, in one shot I have all the output data I want.

outputdata=Table[sol={g1,g2,g3,g4,y} /. FindRoot[ A1[[i]]*(y[[i]])^(k1) + B1[[i]]*(y[[i]])^(k2)==c[[i]] , {y, sol[[5]]+10, sol[[5]], sol[[5]]+20}], {i,1, 60}]

This process worked for a charm for little while but for a certain parameter space (by parameter I don’t mean the xix_i I used earlier but a host of γ\gammas andβ\betas in my equations that I have suppressed so far), it has started giving me errors, 1/0 infinity type stuff, because of some assignment issues. Is there a clean and correct/good way to do this? I wanna be able to import a ton of data, Table my findroot to compute a whole bunch of data and Export it real fast. Please please please help!

MWE

f = y^(3.1276)*(A1) + y^(-0.5875)*(B1) + (c1)^2;
A1 = x1/y + x2*y + 3*x3;
B1 = x1*x3 + 1/(y*x2);
c1 = x1^3 + x2^5 – x3;
g1 = y^(x1) – x3*x2;
g2 = x1/y;

Imported Data Below:

x1 = {89, 88, 87}
x2 = {0.048334203`, 0.048515211`, 0.048707816`}
x3 = {-19486.2273`, -19742.04035`, -20016.22863`}

When I do this, I can see what the curves look like:

Plot[Table[y^(3.1276)*(A1[[i]]) + y^(-0.5875)*(B1[[i]]) + (c1[[i]])^2, {i, 1,
3}], {y, 150, 180}]

This is what I want to do. To be able to Table a whole bunch of output in one shot:

dataoutput=Table[{g1[[i]], g2[[i]], y} /. FindRoot[y^(3.1276)*(A1[[i]]) +
y^(-0.5875)*(B1[[i]]) + (c1[[i]])^2, {y, 165, 150, 170}], {i, 1,
3}]

This is the result that I get, with some errors about accuracygoal and precisiongoal.

{{1623.03, -6.88842*10^19, 167.181}, {1530.37, -5.38632*10^19,
163.049}, {1431., -4.18019*10^19, 158.952}}

My MWE is working, just like my actual problem worked for a certain parameter space but now I am running into trouble. Is there a good way to Import a ton of Data, Table my FindRoot and generate a ton of output using Table and then Export my results? Thanks.

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

  

 

Please including a minimal example in Mathematica code. This is too vague to follow and all the “action” is inside the mystery file that is being imported.
– R. M.♦
Sep 9 ’12 at 20:25

  

 

Thanks RM. I will post a MWE. There is nothing mysterious in the imported file but data. Just 60 values each for x1,x2,…,x5x_1,x_2,…,x_5
– Amatya
Sep 9 ’12 at 20:37

  

 

You can reduce the complexity till you can settle on something small enough that demonstrates your problem. Like 3-4 values instead of 60, etc. (unless if the number 60 is the issue)…
– R. M.♦
Sep 9 ’12 at 20:45

  

 

MWE included with 3 values of Imported data.
– Amatya
Sep 9 ’12 at 21:19

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

1 Answer
1

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

Error Free:

You have almost solved your problem and I will try to help you streamline the whole process involving “Tons of data”! Lets first figure out the RHS of your equation in a symbolic form.

f = y^(3.1276)*(A1) + y^(-0.5875)*(B1) + (c1)^2;
A1 = x1/y + x2*y + 3*x3;r
B1 = x1*x3 + 1/(y*x2);
c1 = x1^3 + x2^5 – x3;
g1 = y^(x1) – x3*x2;
g2 = x1/y;
f // FullSimplify

(x1^3 + x2^5 – x3)^2 + (1 + x1*x2*x3*y)/(x2*y^1.5875) + y^2.1276*(x1 + y*(3*x3 + x2*y))

Now we set the parameter values.

x1 = {89, 88, 87};
x2 = {0.048334203`, 0.048515211`, 0.048707816`};
x3 = {-19486.2273`, -19742.04035`, -20016.22863`};
rhs = Table[y^(3.1276)*(A1[[i]]) + y^(-0.5875)*(B1[[i]]) + (c1[[i]])^2, {i, 1,3}];
scale = rhs /. y -> 150

{1.50954*10^11, 1.12911*10^11, 7.63367*10^10}

Now the trick to make FindRoot free of error message. We know for a an equation f(x)=0f(x)=0 one can scale the RHS by a constant say α∈R\alpha \in \mathbb{R} to get 1αf(x)=˜f(x)=0\frac{1}{\alpha}f(x)=\tilde{f}(x)=0.We will scale each of your three equations by the above scaling factor respectively and then call FindRoot.

Plot[Evaluate[rhs/scale], {y, 150, 180}, Frame -> True,AxesStyle -> Dashed ]

MMA does not spill any error message!

eqs = (# == 0) & /@ (rhs/scale);
res = FindRoot[#, {y, 150}] & /@ eqs

{{y -> 167.181}, {y -> 163.049}, {y -> 158.952}}

However I was not able to reproduce your end result involving FindRoot. Notice that first coordinates are extremely big numbers due to the definition of g1.

dataoutput = Table[{g1[[i]], g2[[i]], y} /. res[[i]], {i, 1, 3}]

{{7.30459*10^197, 0.532358, 167.181}, {4.83167*10^194, 0.539714,
163.049}, {3.23723*10^191, 0.547335, 158.952}}

Export and large data set:

To mimic your real data here I form a larger sample data of size 1000010000. I assume that each parameter x1,x2,x3 follows a NormalDistribution with certain mean (here mean of x1,x2,x3 you supplied) and some arbitrary variance. Then I draw random sample from those distributions.

dataX1 = Floor@RandomVariate[NormalDistribution[Mean[x1], 3], 10^4];
dataX2 = RandomVariate[NormalDistribution[Mean[x2], .5], 10^4];
dataX3 = RandomVariate[NormalDistribution[Mean[x3], 16], 10^4];
TonsOfData = Transpose@{dataX1, dataX2, dataX3};

Now we form a function that does all the calculations and Apply it on TonsOfData. It took just arounf 7.5 seconds in my laptop.

With[{x1 = #[[1]], x2 = #[[2]], x3 = #[[3]]},
Rhs = (x1^3 + x2^5 – x3)^2 + 1 + x1 x2 x3 y/(x2 y^1.5875`) +
y^2.1276` (x1 + y (3 x3 + x2 y));
Scaler = Rhs /. y -> 150;
sol = FindRoot[Rhs/Scaler, {y, 150}];
{y^(x1) – x3*x2, x1/y, y} /. sol
] & /@ TonsOfData; // AbsoluteTiming

{7.4364254, Null}

See how the roots are found gradually by FindRoot using Monitor for a data set of 10001000.

Now we simply export the result into a file called data.csv. Read the documentation to check which other data formats are supported by MMA.

Export[“C:\\Users\\MMA\\Desktop\\data.csv”, TonsOfResult, “Data”];

BR

  

 

PlatoManiac, can I have your babies?
– Amatya
Sep 10 ’12 at 19:53

  

 

Also, if it’s not too much trouble, can you please write in the comments the code that generated the gif with the solutions. The scaling is a great tip, thanks. I’ll try and faithfully reproduce your code for my stuff. Thanks a lot!!!!!!
– Amatya
Sep 10 ’12 at 20:21

  

 

can you please tell me how to extend this command eqs = (# == 0) & /@ (rhs/scale); to the case where rhs/scale is a 2 dimensional list. Thanks. If rhs/scale = {{y,y1},{z,z1}} then right now I am getting {{y,y1}==0,{z,z1}==0} but I would like to get {{y==0,y1==0},{z==0,z1==0}}. Thanks.
– Amatya
Oct 2 ’12 at 21:07