1
$\begingroup$

I'm trying to do a Newthon Raphson with 1st order continuation to solve a problem with 4 variables (fix 2 and obtain the results of 2 (x and Chi).

GaussEliminationwithPivoting[AA_, bb_] := 
  Module[{A = AA, b = bb, n, maxIndex, maxValue, factor, x}, 
   n = Length[b];
   For[i = 1, i <= n, i++, maxIndex = i;
    maxValue = Abs[A[[i, i]]];
    For[j = i + 1, j <= n, j++, 
     If[Abs[A[[j, i]]] > maxValue, maxIndex = j;
       maxValue = Abs[A[[j, i]]];];];
    If[maxIndex != i, A[[{i, maxIndex}]] = A[[{maxIndex, i}]];
     b[[{i, maxIndex}]] = b[[{maxIndex, i}]];];
    For[j = i + 1, j <= n, j++, factor = A[[j, i]]/A[[i, i]];
     A[[j]] -= factor*A[[i]];
     b[[j]] -= factor*b[[i]];];];
   x = Table[0, {n}];
   For[i = n, i >= 1, i--, x[[i]] = b[[i]];
    For[j = i + 1, j <= n, j++, x[[i]] -= A[[i, j]]*x[[j]];];
    x[[i]] /= A[[i, i]];];
   {x}];

(*Functions f1 e f2*)
f1[x_, y_, chi_, v_] := 
  Log[x] - Log[
    y] + (1 - 1/v) (1 - x) - (1 - 1/v) (1 - y) + ((1 - x)^2)/
    chi - ((1 - y)^2)/chi;
f2[x_, y_, chi_, v_] := 
  Log[1 - x] - 
   Log[1 - y] + (1 - v) x - (1 - v) y + (v x^2)/chi - (v y^2)/chi;

(*System of equation*)
equations[x_, chi_, v_, y_] := {f1[x, y, chi, v], f2[x, y, chi, v]};

(*Jacobian*)
jacobian[x_, chi_, v_, y_] := D[equations[x, chi, v, y], {{x, chi}}];

(*Newton Raphson method*)
newtonRaphson[initGuess_List, v_, y_, tol_: 1*^-6, maxIter_: 500] := 
 Module[{X = initGuess, F, J, deltaX, iter = 0}, 
  While[iter < maxIter, F = equations[X[[1]], X[[2]], v, y];
   J = jacobian[X[[1]], X[[2]], v, y];
   deltaX = GaussEliminationwithPivoting[J, -F];
   X = X + deltaX; v = v + deltaX[[-1]];
   If[Norm[deltaX] < tol, Break[]];
   iter++;];
  X]

(*Initial Values*)
v = 1;
initialGuesses = {};

(*Find the initial guess*)
Do[x0 = 1 - y;
 ratio = (1 - x0)/x0;
 chi0 = (1 - (2 x0))/Log[ratio];
 {x, chi} = newtonRaphson[{x0, chi0}, v, y];
 If[0 < x < 1, AppendTo[initialGuesses, {y, x, chi}]], {y, 0.001, 
  0.999, 0.01}]

(*Values of v*)
vVals = Range[50, 500, 50];

(*Find the solution for differents values of v*)
Do[chiVals = {};
 xVals = {};
 yVals = {};
 For[{y, x, chi} = #, {x, chi} = newtonRaphson[{x, chi}, v, y];
    If[0 < x < 1 && x != y, AppendTo[yVals, y]; AppendTo[xVals, x]; 
     AppendTo[chiVals, chi]]] & /@ initialGuesses;
 If[v == 500, 
  ListPlot[Transpose[{yVals, chiVals}], PlotStyle -> Blue, 
   Joined -> True, 
   PlotLabel -> "Plot of chi versus y for v=" <> ToString[v], 
   AxesLabel -> {"y", "chi"}, GridLines -> Automatic]], {v, vVals}]

EDIT:

(*Gauss elimination fuction will be necessary in the code to solve \
the linear system J \[CenterDot] dx = -F *)
GaussEliminationwithPivoting[AA_, bb_] := 
  Module[{A = AA, b = bb, n, maxIndex, maxValue, factor, x}, 
   n = Length[b];
   For[i = 1, i <= n, i++, maxIndex = i;
    maxValue = Abs[A[[i, i]]];
    For[j = i + 1, j <= n, j++, 
     If[Abs[A[[j, i]]] > maxValue, maxIndex = j;
       maxValue = Abs[A[[j, i]]];];];
    If[maxIndex != i, A[[{i, maxIndex}]] = A[[{maxIndex, i}]];
     b[[{i, maxIndex}]] = b[[{maxIndex, i}]];];
    For[j = i + 1, j <= n, j++, factor = A[[j, i]]/A[[i, i]];
     A[[j]] -= factor*A[[i]];
     b[[j]] -= factor*b[[i]];];];
   x = Table[0, {n}];
   For[i = n, i >= 1, i--, x[[i]] = b[[i]];
    For[j = i + 1, j <= n, j++, x[[i]] -= A[[i, j]]*x[[j]];];
    x[[i]] /= A[[i, i]];];
   {x}];


(*Functions f1 and f2*)
f1[x_, y_, chi_, v_] := 
  Log[x] - Log[
    y] + (1 - 1/v) (1 - x) - (1 - 1/v) (1 - y) + ((1 - x)^2)/
    chi - ((1 - y)^2)/chi;
f2[x_, y_, chi_, v_] := 
  Log[1 - x] - 
   Log[1 - y] + (1 - v) x - (1 - v) y + (v x^2)/chi - (v y^2)/chi;

(*Equations*)
equations[x_, chi_, v_, y_] := {f1[x, y, chi, v], f2[x, y, chi, v]};


(*Jacobian*)
jacobian[x_, chi_, v_, y_] := D[equations[x, chi, v, y], {{x, chi}}];

(*Newton-Raphson Method*)
newtonRaphson[initGuess_List, v_, y_, tol_: 1*^-6, maxIter_: 500] := 
 Module[{X = initGuess, F, J, deltaX, iter = 0}, 
  While[iter < maxIter, F = equations[X[[1]], X[[2]], v, y];
   J = jacobian[X[[1]], X[[2]], v, y];
   deltaX = GaussEliminationwithPivoting[J, -F];
   X = X + deltaX;
   If[Norm[deltaX] < tol, Break[]];
   iter++;];
  X]

(*Initial Guess*)
v = 1;
initialGuesses = {};

(*Find the initial guess*)
Do[x0 = 1 - y;
 ratio = (1 - x0)/x0;
 chi0 = (1 - (2 x0))/Log[ratio];
 {x, chi} = newtonRaphson[{x0, chi0}, v, y];
 If[0 < x < 1, AppendTo[initialGuesses, {y, x, chi}]], {y, 0.001, 
  0.999, 0.01}]

(*Valores de v*)
vVals = Range[50, 500, 50];

(*Find solutions for different values of v*)
Do[chiVals = {};
 xVals = {};
 yVals = {};
 For[{y, x, chi} = #, {x, chi} = newtonRaphson[{x, chi}, v, y];
    If[0 < x < 1 && x != y, AppendTo[yVals, y]; AppendTo[xVals, x]; 
     AppendTo[chiVals, chi]]] & /@ initialGuesses;
 If[v == 500, 
  ListPlot[Transpose[{yVals, chiVals}], PlotStyle -> Blue, 
   Joined -> True, 
   PlotLabel -> "Plot of chi versus y for v=" <> ToString[v], 
   AxesLabel -> {"y", "chi"}, GridLines -> Automatic]], {v, vVals}]

And I'm obtained this errors and the code dont compile.

During evaluation of In[12]:= General::ivar: 0.999` is not a valid variable.

During evaluation of In[12]:= Thread::tdlen: Objects of unequal length in {{0.999,0.144496}}-{0.999,-0.0722481} cannot be combined.

During evaluation of In[12]:= Part::partd: Part specification {0,0}[[2,2]] is longer than depth of object.

During evaluation of In[12]:= Part::partd: Part specification {0,0}[[2,2]] is longer than depth of object.

During evaluation of In[12]:= Part::partd: Part specification {0,0}[[1,2]] is longer than depth of object.

During evaluation of In[12]:= General::stop: Further output of Part::partd will be suppressed during this calculation.

During evaluation of In[12]:= Thread::tdlen: Objects of unequal length in {0.999,0.144496}+{{{(-1.7763610^-15-0.999 Part[<<3>>] Power[<<2>>])/{<<2>>}[[1,1]],(-1.7763610^-15-0.144496 Part[<<3>>] Power[<<2>>])/{<<2>>}[[1,1]]},{0.999/{<<2>>}[[2,2]],0.144496/{<<2>>}[[2,2]]}}} cannot be combined.

During evaluation of In[12]:= Norm::nvm: The first Norm argument should be a scalar, vector, or matrix.

During evaluation of In[12]:= Thread::tdlen: Objects of unequal length in {{{(1-Power[<<2>>] Plus[<<2>>])^2,(1-Power[<<2>>] Plus[<<2>>])^2},{(1-0.999 Power[<<2>>])^2,(1-0.144496 Power[<<2>>])^2}}} {1.001,6.9206} cannot be combined.

During evaluation of In[12]:= General::stop: Further output of Thread::tdlen will be suppressed during this calculation.

During evaluation of In[12]:= General::ivar: {{{(-1.7763610^-15-(0.999 {<<2>>}[[1,2]])/Part[<<3>>])/{0,0}[[1,1]],(-1.7763610^-15-(0.144496 {<<2>>}[[1,2]])/Part[<<3>>])/{0,0}[[1,1]]},{0.999/{0,0}[[2,2]],0.144496/{0,0}[[2,2]]}}} is not a valid variable.

During evaluation of In[12]:= Norm::nvm: The first Norm argument should be a scalar, vector, or matrix.

During evaluation of In[12]:= General::ivar: {{{{{{Power[<<2>>] Plus[<<6>>],Power[<<2>>] Plus[<<6>>]},{Power[<<2>>] Plus[<<6>>],Power[<<2>>] Plus[<<6>>]}}},{(-6.90776+{{<<2>>}}+{{<<2>>}}+{0.999,<<18>>}-{<<1>>} {<<2>>}-Plus[<<6>>] Part[<<3>>] Power[<<2>>])/{<<2>>}[[1,1]],(-6.90776+<<7>>)/{<<2>>}[[1,1]]}},{{{{Plus[<<6>>] Power[<<2>>],Plus[<<6>>] Power[<<2>>]},{Plus[<<6>>] Power[<<2>>],Plus[<<1>>] <<1>>}}},<<1>>}}} is not a valid variable.

During evaluation of In[12]:= General::stop: Further output of General::ivar will be suppressed during this calculation.

During evaluation of In[12]:= Part::partw: Part 2 of {{{{{{{{<<2>>},{<<2>>}}},{Times[<<4>>]+Times[<<2>>],Times[<<4>>]+Times[<<2>>]}},{{{{<<2>>},{<<2>>}}},{Times[<<4>>]+Times[<<2>>],Times[<<4>>]+Times[<<2>>]}}}},{{{-{<<1>>} {<<2>>} {<<2>>}+Power[<<2>>] Plus[<<2>>],-{<<1>>} {<<2>>} {<<2>>}+Power[<<2>>] Plus[<<2>>]},{-{<<1>>} {<<2>>} {<<2>>}+0.999 Power[<<2>>],-{<<1>>} {<<2>>} {<<2>>}+0.144496 Power[<<2>>]}}}}} does not exist.

During evaluation of In[12]:= Part::partw: Part 2 of {{{{{{{{<<2>>},{<<2>>}}},{Times[<<4>>]+Times[<<2>>],Times[<<4>>]+Times[<<2>>]}},{{{{<<2>>},{<<2>>}}},{Times[<<4>>]+Times[<<2>>],Times[<<4>>]+Times[<<2>>]}}}},{{{-{<<1>>} {<<2>>} {<<2>>}+Power[<<2>>] Plus[<<2>>],-{<<1>>} {<<2>>} {<<2>>}+Power[<<2>>] Plus[<<2>>]},{-{<<1>>} {<<2>>} {<<2>>}+0.999 Power[<<2>>],-{<<1>>} {<<2>>} {<<2>>}+0.144496 Power[<<2>>]}}}}} does not exist.

During evaluation of In[12]:= Norm::nvm: The first Norm argument should be a scalar, vector, or matrix.

During evaluation of In[12]:= General::stop: Further output of Norm::nvm will be suppressed during this calculation.

Out[20]= $Aborted

Someone could help me please ?

$\endgroup$
2
  • 1
    $\begingroup$ When computing you initial guess, you call newtonRaphson with {x0,chi0} = {.999,.144..}. These get assigned to X, and then to equations[] and jacobian. jacobian should probably not have list braces for the last arg, but a more serious problem is that you are plugging numbers in for x,chi and trying to compute derivative wrt .999,.144. These need to be variables/symbols in order to compute derivatives. I don't understand your algorithm, so you'll have to decide how to fix this. A couple suggestions. Mathematica is an interpreted (not compiled) language, so you can test your individual functions $\endgroup$
    – user87932
    Oct 10, 2023 at 2:49
  • 1
    $\begingroup$ ..interactively. Since you have set tolerance and maxIter, make them small while debugging. For newRaph, I set {tol,maxIter} = {1,2} to avoid hanging up in the loop, and set the step for y in your initGuess Do loop to 1. Finally, the debugger in the evaluation menu by default breaks on messages/assertions (which you can add). Try launching this when you run; it halts on the first message and shows you an evaluation stack which may give you helpful information. If you can fix the problem I mentioned above, I can take another look, but you'll have to decide what the code should be doing. $\endgroup$
    – user87932
    Oct 10, 2023 at 2:53

1 Answer 1

2
$\begingroup$

A description of the Newton-Raphson method for the scalar and k variable, k equation case may be found here https://en.wikipedia.org/wiki/Newton%27s_method

For a scalar function, the roots of f(x) may be found by iterating x[n+1] = x[n] f[x[n]]/f'[x[n]].

The matrix form of the Newton-Raphson method for a vector x can be computed by solving the system of linear equations given by

$J_F(x_n) (x_{n+1}- x_n) == -F(x_n)$

where x is a vector, F is a vector valued function, and $J _F(x)$ is the Jacobian of F. Examining the newtonRaphson function in the post, F is computed, then a Jacobian, deltas for x and v are computed by Gaussian elimination, then the deltas are used to update the x,v variables. The rest of the code controls the iteration. That makes sense in view of the equation we're iterating over.

It should be possible to simplify the code by taking advantage of the linear algebra functions provided by Mathematica. There seem to be problems with the code for Gaussian elimination. In an earlier comment, I mentioned that the various functions can be tested individually, which is what I did. I tried manually entering an A and b matrix into the Gaussian elimination code and got an output with Null and numerical values which don't seem to be correct.

enter image description here

Mathematica's built-in LinearSolve does produce the correct solution, as can easily be verified by substituting the result back in and evaluating A.x ==b

enter image description here

So I'm going to post a solution which takes advantage of LinearSolve . Mathematica utilizes widely used standard libraries such as LAPACK, etc. for linear algebra, so LinearSolve is a reliable substitute.

EDIT: Due to an error on my part, the code in newtonRaphson wasn't running; removing the Break[] caused it to immediately exit since the test compared a symbol to a number. After fixing this, I found there were a number of problems with the Jacobian calculation, which I corrected.

I split the argument lists up into parameters (set once for each run and never changed) and variables which are initialized and updated while the newtonRaphson algorithm runs. Otherwise, it's hard for someone unfamiliar with this code to understand which is which without reading it carefully.

Edit: I got the order swapped around: y,v instead of v,y. This cured the singular matrix problem I saw earlier.

(*Functions f1 and f2. Treat v,y as parameters to be set to constant \
values and x,chi as variables. *)
   f1[v_,y_][x_, chi_] := 
  Log[x] - 
   Log[y] + (1 - 1/v) (1 - x) - (1 - 1/v) (1 - y) + ((1 - x)^2)/
    chi - ((1 - y)^2)/chi;
f2[v_,y_][x_, chi_] := 
  Log[1 - x] - 
   Log[1 - y] + (1 - v) x - (1 - v) y + (v x^2)/chi - (v y^2)/chi;

(*Equations*)
equations[param1_, param2_][var1_, 
   var2_] := {f1[param1, param2][var1, var2], 
   f2[param1, param2][var1, var2]};

The next problem is the handling of the Jacobian. What we want to do is to first compute a symbolic Jacobian from the equations, then use that for the rhs of a function which replaces the symbolic arguments with numbers to get a numerical matrix.

EDIT: I didn't check the original equation for the Jacobian carefully. You need to use an Outer product to compute it. I use dummy variables to compute it symbolically, then substitute the numerical values into the symbolic form.

(*Jacobian*)
numericJacobian[param1_?NumericQ, param2_?NumericQ][varVal1_?NumericQ,
    varVal2_?NumericQ] :=
    Module[{dummy1, dummy2},
        Outer[D, 
     equations[param1, param2][dummy1, dummy2], {dummy1, 
      dummy2}] /. {dummy1 -> varVal1, dummy2 -> varVal2}
    ];

Replace Gaussian elimination with LinearSolve and pass numeric parameters into the symbolic Jacobian to get a numerical matrix. I fixed the loop exit code which I'd broken, and modified the code to handle the new function interfaces.

EDIT: I found that the messages about singular matrices were caused by x going negative, then complex during iteration, and added some exception handling/validity indicators to deal with that. I'll comment on this more below, and add some plots to show what's going on.

(*Newton-Raphson Method*)
newtonRaphson[param1_,param2_][varVal1_,varVal2_,tol_:1*^-6,maxIter_:500] :=
    Module[{var1,var2,F,J,delta,iter=0,exit=False,valid="Valid"},
        {var1,var2}={varVal1,varVal2};
        While[!exit,
            F = equations[param1,param2][var1,var2];
            Enclose[
                J = ConfirmQuiet[numericJacobian[param1,param2][var1,var2], All, valid="Failed Jacobian calc"];
                delta = ConfirmQuiet[LinearSolve[J,-F], All, valid="Failed LinearSolve"];
            ];
            {var1,var2} += delta;
            iter++;
            exit = iter >= maxIter || Norm[delta] <= tol; ];
        If[valid==="Valid", {var1,var2,valid}, {varVal1,varVal2,valid}]
    ];

EDIT: The OP requested help in generating plots, so I made some further revisions to the code. One problem is in how "v" was handled. Although a list of v's was generated, it wasn't used. Only v=1 was computed, but the plot was wrapped in test for v==500. Here's the updated version, with an explanation of the changes I made.

I replaced the code for "v" with vVals set to all the values used originally. I.e. 1, along with 50 to 500 in steps of 50. The Do loop had a single iterator over y; I added a second iterator over vVals. The original version ran newtonRaphson twice, with the values for {x,chi} from the first pass. I removed the second call. I renamed "initialGuesses" to "results", since some of the initial values get overwritten by the output of the iteration. Since y,v are needed for plotting, I moved them outside the loop, but iterate over their contents.

EDIT: I removed the test for appending the output to results and added a validity indicator. The consequence of having x stray outside some range is likely to be a singular matrix, and the Confirm/Enclose exception handler should trap that and set the validity indicator which I added to the output of newtonRapshson. I'm going to save everything so you have an idea what inputs cause failures. The data can be filtered afterwards for plots, etc.

(* Run the algorithm. Use y,v as parameters, x,chi as variables which \
are set to initial guesses, then updated.*)
results = {}; 
vVals = Join[{1}, Range[50, 500, 50]]; yVals = Range[.001, .999, .01];
Do[ x0 = 1 - y;
    ratio = (1 - x0)/x0;
    chi0 = (1 - (2 x0))/Log[ratio];
    {x,chi,valid} = newtonRaphson[v, y][x0, chi0, 1*^-6, 10];
    AppendTo[results, {y, v, x, chi, valid}];
  , {y, yVals}, {v, vVals}];

For the plots, I removed the test for v=500 and did plots for all v values. I replaced ListPlot with ListLinePlot since the points were being joined to form a line anyway. I generated the plots using Table to iterate over vVals. I used a Case statement to select data matching specific values of v (the third item in each list), and if a match occurs, transform the matching sublist to {y,chi}. Since I get a number of plots for different v values, I Partition the list of plots into rows of 3; UpTo allows the last sublist to be shorter than 3. I wrap the output in GraphicsGrid. The grid may default to a small size and be hard to read. Just click on the plot and drag the lower right corner down to expand it.

Originally, these plots were of {y,chi} for different v. I thought that odd, because y is treated as a fixed parameter value. It seemed like the plots should show {x,chi} for different {v,y} choices, which is what I did here. I used StringTemplate to set up the label, and reduced the number of y,v values since the result was too big to fit on my screen if I allowed all values to be plotted.

EDIT: Although I put validity flags in to tag situations where error messages are generated, it turns out that you can still have cases which produce no messages but do yield complex results. I added checks to the pattern for x,chi to filter those events out. With these changes, you now see variations in x,chi as you vary the other parameters. The cases where you exit quickly are simply due to the fact that the initial guesses are quite close to the solution.

(*Generate plots of chi vs. x for different y,v values.Partition the \
plots so they form rows with three plots per row.*)
label = StringTemplate["Plot of chi versus x for v= `` and y= ``"];
GraphicsGrid[
    Partition[Flatten[Table[
            ListLinePlot[
     Cases[results, {_, _?(# == v &), x_, chi_} :> {x, chi}],
                PlotStyle -> Blue,
                PlotLabel -> label[v, y],
                AxesLabel -> {"x", "chi"},
                GridLines -> Automatic],
        {y, {.001, .421, .991}}, {v, {1, 50, 250}}]], UpTo[3]]]

enter image description here

This code can be used to filter out results with valid and non-complex entries. I limited the output to 20 entries. To get individual v,y,x,chi,validity values, use Transpose on the output of the Case statement and assign to {yVal,vVal, etc.}.

Take[ Cases[results, {_, _, x_Real, chi_Real, "Valid"} ], 20]

{{0.001, 1, 0.999, 0.144496, "Valid"}, {0.001, 250, 0.00110143, 
  0.00210108, "Valid"}, {0.001, 300, 0.00110554, 0.00210513, 
  "Valid"}, {0.001, 350, 0.00108866, 0.00208821, "Valid"}, {0.001, 
  400, 0.000989439, 0.00198849, "Valid"}, {0.011, 1, 0.989, 0.217391, 
  "Valid"}, {0.011, 50, 0.011955, 0.0229472, "Valid"}, {0.011, 100, 
  0.0120238, 0.0230181, "Valid"}, {0.011, 150, 0.012043, 0.023038, 
  "Valid"}, {0.011, 200, 0.012052, 0.0230474, "Valid"}, {0.011, 250, 
  0.0120573, 0.0230529, "Valid"}, {0.011, 300, 0.0120607, 0.0230565, 
  "Valid"}, {0.011, 350, 0.0120631, 0.023059, "Valid"}, {0.011, 400, 
  0.0120649, 0.0230609, "Valid"}, {0.011, 450, 0.0120663, 0.0230624, 
  "Valid"}, {0.011, 500, 0.0120675, 0.0230635, "Valid"}, {0.021, 1, 
  0.979, 0.249349, "Valid"}, {0.021, 50, 0.0227425, 0.0437187, 
  "Valid"}, {0.021, 100, 0.0227993, 0.0437848, "Valid"}, {0.021, 150, 
  0.0228168, 0.0438055, "Valid"}}

EDIT: Some discussion about the singularities. After making another change to newtonRapshson to fix one of my errors, I began getting error messages again. The three red dots to the left of the message open a short menu when you click on them, with one option being show the evaluation stack. I was able to use this to get the exact parameters used when the error occurred: v=100,y= .001,x=.999,chi- 0.14449622608365997. I then ran newtonRaphson stand alone, with some additional outputs added to aid in debugging. I found that the problem occurred quickly, so I extracted the three lines to compute the equations, the Jacobian, and run LinearSolve, and executed them one by one. What happens is that the first pass produces a delta which results in x becoming negative. Repeating with the adjusted x,chi values results in complex numbers being generated, and it goes downhill from there. I looked at the equations, substituted the parameters, computed the Jacobian, examining and checking each stage, and everything is in order. So this is the behavior of the system, not another bug in the code. I plotted your equations with the values for v,y plugged in:

Plot3D\[5.918745278982136` - 0.998001`/chi + 0.99` (1 - x) + (1 - x)^2/
  chi + Log\[x\], {x, 0, 1},][5]][5] {chi, 0, 1}]

Plot3D[0.10000050033358354` - 0.00009999999999999999`/chi - 99.` x + (
  100.` x^2)/chi + Log[1 - x], {x, 0, 1}, {chi, 0, 1}]

GraphicsColumn[{%%, %}]

enter image description here

You have problems near the edges and a relatively flat surface elsewhere for both equations. The behavior near the edges is the cause of your singular matrices.

There's nothing I can do about the behavior of the equations, but I did add exception handling to the code to trap the messages and do an early exit with some sort of "failure" indication - though you can still generate complex outputs without any messages, which need to be filtered out of the data.

$\endgroup$
2
  • $\begingroup$ Thank you so much for your help, but I tried to plot the graph with this modifications and still having problems (Initial Guess) initialGuesses = {}; (Find the initial guess) Do[x0 = 1 - y; ratio = (1 - x0)/x0; chi0 = (1 - (2 x0))/Log[ratio]; {x, chi} = newtonRaphson[{x0, chi0}, v, y]; If[0 < x < 1, AppendTo[initialGuesses, {y, x, chi}]], {y, 0.01, 0.99, 0.01}] (v values) vVals = Range[50, 500, 50]; $\endgroup$ Oct 12, 2023 at 16:50
  • $\begingroup$ chiResults = {}; (find the v) Do[chiVals = {}; ({y, x, chi} = #; {x, chi} = newtonRaphson[{x, chi}, v, y]; If[0 < x < 1 && x != y, AppendTo[chiVals, chi]]) & /@ initialGuesses; AppendTo[chiResults, chiVals], {v, vVals}] (Create the plot of chi versus y for all v values) ListLinePlot[Transpose[chiResults], PlotRange -> All, PlotLabel -> "Plot of chi versus y for all v values", AxesLabel -> {"y", "chi"}, GridLines -> Automatic, PlotLegends -> vVals] $\endgroup$ Oct 12, 2023 at 16:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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