Examples don't match codebase

[mludlum]

I'd like to use your library, but the the code in this hub does not support the examples you give. Do you have an update to your code that has the IFunction interface and supporting classes?

[EduardBargues]

Hi, mludlum. Thank you for opening an issue. The code of this repository has changed quite a bit since the last time I updated. Let me put the last version and check if it is still suitable for you. I'd say that in 48 hours I'll be able to update this repository. Thank you for your patience!

[mludlum]

That would be great. I like the framework in your examples versus the test case code in the repo. In the meantime, I'm experimenting with this I found on Rosetta code: https://rosettacode.org/wiki/Multidimensional_Newton-Raphson_method#C.23

I'd prefer to use MathNet for the data structures like you do.

Are you considering creating a pull request to add your code to MathNet? It seems that their solvers are only for linear systems.

[EduardBargues]

Uau! Those are big words :)! I have never thought about creating a PR to MathNet. I think my code would need some extra "style" work to fit their requirements. MathNet does not provide a solver like the one in this repo. I'll think about creating a pr to their repo.

[EduardBargues]

Hi mludlum ! I just updated the code and the readme of the repo. Let me know if you miss anything or have some other questions. Have a great day!

[mludlum]

Thanks for the update. Please add Program.cs to the Test project. https://github.com/EduardBargues/NonLinearEquationsSolver/tree/master/Test.NonLinearEquationsSolver

[mludlum]

I'm not getting good results, so I must be doing something wrong. I'd appreciate any suggestions. I don't know if newton-raphson is the right approach to my problem. When I tried this using rosetta code (https://rosettacode.org/wiki/Multidimensional_Newton-Raphson_method#C.23) it never converged so I thought I'd try your code. However, I don't understand the terms force, energy, load, etc. Can you suggest values?

I'm trying to solve the problem of getting calibration offsets for 3 angle sensors on an excavator that has 3 arms. I'm using forward kinematics to compare the measured reach of the excavator (the distance projected on the ground from the start of arm1 to the end of arm3). My function and derivative below contains sines and cosines. L[j] is the length of each of the 3 arms. A[j,k] is the measured angle and X[k] is the measured reach where j=0..2 for each of the arms and k=0..2 for each of the three sets of measurements I'm taking. I don't have a good initial guess of the 3 offsets so I'm using zeros. My units are radians, so I'd like my result to be within 0.001 accuracy.

`
Vector Function(Vector offsets) { Vector result = new DenseVector(3); for (int j = 0; j < 3; j++) { result[j] = L[0] * Math.Cos(Math.Abs(A[0, j] + offsets[0]))

• L[1] * Math.Cos(Math.Abs(A[1, j] + offsets[1]))
  • L[2] Math.Cos(Math.Abs(A[2, j] + offsets[2])) - X[j]; } return result; } Matrix Stiffness(Vector offsets) { Matrix matrix = new DenseMatrix(3, 3); for (int j = 0; j < 3; j++) { for (int k = 0; k < 3; k++) { matrix[j, k] = -L[j] Math.Sin(Math.Abs(A[j, k] + offsets[j])); } } return matrix; } DenseVector force = new DenseVector(3) { [0] = 1, [1] = 1, [2] = 1 }; double inc = 0.0001; Solver solver = Solver.Builder .Solve(3, Function, Stiffness) .Under(force) .UsingStandardNewtonRaphsonScheme(inc) //.WithInitialConditions(0.001, DenseVector.Create(3, 0), DenseVector.Create(3, 1)) .Build(); List states = solver.Broadcast().TakeWhile(x => x.Lambda <= 0.0001).ToList();

`

[EduardBargues]

Hi again! About the test project I just remove it. It was something that I was not using. The good tests are inside the main project in the classes TestUtils, TestSolverResults and TestSolverBuilder. About your problem, I'll take a look and let you know.

[EduardBargues]

Hi! I read your problem and it is going to need more explanation. I'm going to close this issue because it changes the subject quite fast and open a new one with your question. See you there! Okay, I though it was easier to open new issues but I don't know how to tag you in another one. So I'll answer here. There we go! :)

• First of all, this repository is not completely mature. It served me well for my thesis but obviously I did not tests all the possible scenarios. Is still under development and I would like to do more tests (tones of them).
• It is a good question to ask if the Newton-Raphson is the best method for your system. The answer is yes but you may need a good initial aproximation of try a more powerful pseudo Newton-Raphson method such as Arc-Length.
• About the terms you asked for. Esentially we are solving a system like R(u)=lambda*f, where R is the "reaction" of our system, f is the reference external load, lambda is our loadFactor (usually goes from 0 to 1) and u is the displacements (degrees of freedom of our system. When you try to solve a system of this kind, you want the iterative method to converge in terms of three tolerance: displacements (stop when the solution obtained does not change much), force (stop when your Reaction counteracts your external force), and energy (stop when work done by your reaction is equal to the work done by the external forces along the equilibrium path. In all 3 tolerances, you should be okay setting it's value to 1e-3 (I'm adding a pdf of a chapter of my thesis where you can see the equations and how they work).
• About providing predifined values, there is no straight answer. Sorry 👍 ... Let me get back to it tomorrow and I'll try to solve it and let you know.

[EduardBargues]

This is your code commented with some tips and my proposal. Let me know if it works :) !

Vector Function(Vector offsets)
{
    Vector result = new DenseVector(3);
    for (int j = 0; j < 3; j++)
    {
        result[j] = L[0] * Math.Cos(Math.Abs(A[0, j] + offsets[0]))
            + L[1] * Math.Cos(Math.Abs(A[1, j] + offsets[1]))
            + L[2] * Math.Cos(Math.Abs(A[2, j] + offsets[2])) - X[j];
    }
    return result;
}
Matrix Stiffness(Vector offsets)
{
    Matrix matrix = new DenseMatrix(3, 3);
    for (int j = 0; j < 3; j++)
    {
        for (int k = 0; k < 3; k++)
        {
            matrix[j, k] = -L[j] * Math.Sin(Math.Abs(A[j, k] + offsets[j]));
        }
    }
    return matrix;
}
DenseVector force = new DenseVector(3) { [0] = 1, [1] = 1, [2] = 1 };
double inc = 0.0001; // Here you are setting a lambda increment equal to 1e-4 !! Really small in my opinion.
Solver solver = Solver.Builder
    .Solve(3, Function, Stiffness) // three dof, perfect. The Function and Stiffness seem to be correctly defined.
    .Under(force) // External force equal to this vector -> (1,1,1)
    .UsingStandardNewtonRaphsonScheme(inc) // You are using a standard newton-raphson approach with lambda increment = 1e-4. 
    //This means that the solver will go slowly from an external force (0,0,0) to (1,1,1) incrementing the external by inc*(1,1,1)
    // on each step. On each step, it will launch a NR algorithm to equilibrate the system in this new configuration.
    .Build();
List states = solver.Broadcast().TakeWhile(x => x.Lambda <= 0.0001).ToList(); // Here you are asking the solver to provide
// you with the intermediate equilibrated stats untill it reaches a load factor of 0.0001. But I think that you want to find the
// solution when the system is under the (1,1,1) external load, right?

// If you think your solver can find the solution in one step I recommend you to build a different solver as I provide you in this 
// piece of code:
Solver.Builder
    .Solve(3,Function, Stiffness)
    .Under(force)
    .UsingStandardNewtonRaphsonScheme(1) // Find solution in one single step. the solver will equilibrate just two times the system:
    // one under lambda=0 and one with lambda=1
    .Build();
var solution = solver.Broadcast().TakeUntil(x=>Math.Abs(1-x.Lamda)<=1e-3).Last(); // I'm taking the last solution of all intermediate solutions (TakeUntil is not in linq, use morelinq)

[EduardBargues]

Closing issue since no further answer has been provided :) ...