forum.alglib.net

ALGLIB forum
It is currently Tue Feb 27, 2024 12:14 am

All times are UTC


Forum rules


1. This forum can be used for discussion of both ALGLIB-related and general numerical analysis questions
2. This forum is English-only - postings in other languages will be removed.



Post new topic Reply to topic  [ 11 posts ]  Go to page 1, 2  Next
Author Message
 Post subject: nonlinear equations system
PostPosted: Sat Jun 04, 2011 4:19 pm 
Offline

Joined: Sat Jun 04, 2011 3:58 pm
Posts: 5
I have a system of non-linear equations which has to be solved.
This problem was previously solved by means of mathcad, but now it's considered unconvinient and I have to program an application with cute GUI etc. But, the problem is, I can't find any way to solve the system.
Mathcad claims it uses levenberq-marquardt method to solve my system, so, obviously, I tried to use Alglib's minlm* routines. But the algorithm stops because of "relative step is no more than EpsX"(terminationtype == 2), and the solution provided is not as good as I want it to be (and not as Mathcad's solution is).
So I'm looking for any advice on how to get a proper result (or which other algorithm to use).
Thanks in advance.


Top
 Profile  
 
 Post subject: Re: nonlinear equations system
PostPosted: Sun Jun 05, 2011 7:13 pm 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 900
Can you post your equation here along with starting point? And, if you can, your code which calls ALGLIB. We can try to find the reason for the failure to converge.

I am really interested in solving this problem because I want ALGLIB to be as good as possible.


Top
 Profile  
 
 Post subject: Re: nonlinear equations system
PostPosted: Mon Jun 06, 2011 1:17 pm 
Offline

Joined: Sat Jun 04, 2011 3:58 pm
Posts: 5
the original system is in .png file attached to this reply
Attachment:
File comment: system
system.png
system.png [ 22.22 KiB | Viewed 8352 times ]

unknown variables are L, H, A and R_x. other variables will be taken from the input but now I treat them as constants for testing purpose.
the code is as follows
Code:
using System;
using System.Text;

namespace Project1
{
    class Program
    {
        //init
        static double d = 3;
        static double a = 1000;
        static double F_0 = 0.035;
        static double F_H = 549;
        static double F_d = 9;
        static double c_yd = 1.3;
        static double c_xd = 0.5;
        static double F_sh = 5;
        static double c_ysh = 0.6;
        static double c_xsh = 0.9;
        static double M_g = 1500;
        static double M_z = 400;
        static double hamma_M = 74000;
        static double V = 2.2;
        //
        static double ro = 1030;
        static double v = 1.3e-6;
        static double g = 9.80665;
        static double d_a = d / a;
        static double F_f = F_H * a / d;
        static double G_g = M_g * g * (hamma_M - ro * g) / hamma_M;
        static double G_z = M_z * g * (hamma_M - ro * g) / hamma_M;
        static double R_yd = c_yd * ro * V * V * F_d / 2;
        static double R_xd = c_xd * ro * V * V * F_d / 2;
        static double Q = c_ysh * ro * V * V * F_sh / 2;
        static double R_xsh = c_xsh * ro * V * V * F_sh / 2;
       
        //
        static void f(double[] x, double[] fi, object obj)
        {
            double L = x[0];
            double H = x[1];
            double A = x[2];
            double R_x = x[3];
            fi[0] = R_x - (0.04 * (A + 110.0 * Math.Exp(-6.74 * F_0) * (Math.PI * H * L / (4.0 * F_f) - 0.02)) - 0.09) * (ro * V * V * F_H / 2.0);

            fi[1] = Math.PI * H * L / (4.0 * F_f) - 1.0 + Math.Exp(0.0008 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / ((R_x + R_xsh) * F_0));

            fi[2] = A + 110.0 * Math.Exp(-6.74 * F_0) * (Math.PI * H * L / (4.0 * F_f) - 0.02) - A * (1.0 - Math.Exp(-4.5 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / (R_x + R_xsh)));

            fi[3] = H / L - Math.Exp(1.0 - 2.0 * R_yd / (Q + 2.0 * G_g + G_z));
        }

        static void Main()
        {




            double[] x = { 50, 50, 6, 10000 };

            double[] fi = { 0, 0, 0, 0 };

            double epsg = 0.0000000000001;
            double epsf = 0;
            double epsx = 0;
            int maxits = 0;
            alglib.minlmstate state;
            alglib.minlmreport rep;

            alglib.minlmcreatev(4, x, 0.0001, out state);
            alglib.minlmsetcond(state, epsg, epsf, epsx, maxits);
            alglib.minlmoptimize(state, f, null, null);
            alglib.minlmresults(state, out x, out rep);

            System.Console.WriteLine("{0}", rep.terminationtype); // EXPECTED: 4
           
            f(x, fi, null);
            System.Console.WriteLine(" X : {0}, {1}, {2}, {3}", x[0], x[1], x[2], x[3]);
            System.Console.WriteLine(" f(X) : {0}, {1}, {2}, {3}", fi[0], fi[1], fi[2], fi[3]);
            System.Console.ReadLine();
        }
    }
}



it's almost completely copied from alglib example for minlv*, save for the function and dimension chanched
I sincerely hope, there is any hint on what to do
Anyway, thank you for your responsiveness


Top
 Profile  
 
 Post subject: Re: nonlinear equations system
PostPosted: Tue Jun 07, 2011 10:51 am 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 900
I suggest you to scale your variables with alglib.minlmsetscale() function. Appropriate scale can be [10, 10, 10, 10000]. Unfortunately, I've been too busy yesterday to try it out, but I hope that I will find some time today.


Top
 Profile  
 
 Post subject: Re: nonlinear equations system
PostPosted: Tue Jun 07, 2011 5:40 pm 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 900
I've tested your code. The problem was caused by different scale of the variables, which: a) caused numerical differentiation to lose precision, b) caused incorrect work of the sopping criteria, c) led to incorrect choice of the damping parameter.

Setting scale alglib.minlmsetscale() solved this problem. I also recommend you to use less stringent stopping condition - epsg=1.0E-6 should be enough.


Top
 Profile  
 
 Post subject: Re: nonlinear equations system
PostPosted: Wed Jun 08, 2011 2:49 am 
Offline

Joined: Sat Jun 04, 2011 3:58 pm
Posts: 5
I'm really sorry for being stupid and bothering you, but:
I tried to use minlmsetscale in such a way:
Code:
minlmsetscale(state, new double[]{10,10,10,10000});

right before invoking minlmoptimize(...)
but it didn't help - the solution is still not like the expected one. perhaps, i did sth wrong?
By the way, I'd like to know how to choose an appropriate vector for scaling(at least intuitively)
Thank you for your help


Top
 Profile  
 
 Post subject: Re: nonlinear equations system
PostPosted: Wed Jun 08, 2011 4:39 am 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 900
1. Scale of the variable is its most typical value. You should know from problem formulation that one variable should be about 0.1-10 in magnitude while other can be as large as 100000.

2. Can you post exact (expected) answer? Your code at my computer returns

Quote:
X : 1,48064959286412, 0,805990620146673, 8,63291052919714, 254292,030484443
f(X) : 5,82076609134674E-11, 0,00818196564353979, 1,78589132451279E-05, -4,8129
8789711815E-08


and now I see that x[0] and x[1] are quite non-stable (they wander from 0 to their maximum values), and f[1] is significantly non-zero.

P.S. Don't worry to ask questions, I really appreciate feedback like yours.


Top
 Profile  
 
 Post subject: Re: nonlinear equations system
PostPosted: Wed Jun 08, 2011 5:47 am 
Offline

Joined: Sat Jun 04, 2011 3:58 pm
Posts: 5
I got almost the same result - and f[1]'s not being zero is definitely wrong
thanks to mathcad, the expected answer is (approximately):

Code:
f[0] = 89.84
f[1] = 48.9
f[2] = 4.22
f[3] = 102374


Top
 Profile  
 
 Post subject: Re: nonlinear equations system
PostPosted: Wed Jun 08, 2011 5:31 pm 
Offline

Joined: Sat Jun 04, 2011 3:58 pm
Posts: 5
I've found something, that I consider to be a root of the problem.
I tried to calculate f for the approximate x (given in the message above) and have got an interesting result:

f(x) == 42.52 ; 0.038 ; almost 0 ; almost 0;

hey! f(x)[0] is not 0! but the solution is almost exact! it's not surprising that minlm* jumps over the correct solution and looks for it in wrong direction.
That led me to the idea of scaling f(x) too. The problem is (the new one), I don't know how to do it.
First of all I tried to divide f[i] by left part of corresponding equation, but it didn't help.
The source code now looks like
Code:
using System;
using System.Text;

namespace Project1
{
    class Program
    {
        //init
        static double d = 3;
        static double a = 1000;
        static double F_0 = 0.035;
        static double F_H = 549;
        static double F_d = 9;
        static double c_yd = 1.3;
        static double c_xd = 0.5;
        static double F_sh = 5;
        static double c_ysh = 0.6;
        static double c_xsh = 0.9;
        static double M_g = 1500;
        static double M_z = 400;
        static double hamma_M = 74000;
        static double V = 2.2;
        //
        static double ro = 1030;
        static double v = 1.3e-6;
        static double g = 9.80665;
        static double d_a = d / a;
        static double F_f = F_H * a / d;
        static double G_g = M_g * g * (hamma_M - ro * g) / hamma_M;
        static double G_z = M_z * g * (hamma_M - ro * g) / hamma_M;
        static double R_yd = c_yd * ro * V * V * F_d / 2;
        static double R_xd = c_xd * ro * V * V * F_d / 2;
        static double Q = c_ysh * ro * V * V * F_sh / 2;
        static double R_xsh = c_xsh * ro * V * V * F_sh / 2;
       
        //
        static void f(double[] x, double[] fi, object obj)
        {
            double L = x[0];// *x[0];
            double H = x[1];// * x[1];
            double A = x[2];// * x[2];
            double R_x = x[3];// * x[3];
            fi[0] = (R_x - (0.04 * (A + 110.0 * Math.Exp(-6.74 * F_0) * (Math.PI * H * L / (4.0 * F_f) - 0.02)) - 0.09) * (ro * V * V * F_H / 2.0));

            fi[1] = (Math.PI * H * L / (4.0 * F_f) - 1.0 + Math.Exp(0.0008 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / ((R_x + R_xsh) * F_0)));

            fi[2] = A + 110.0 * Math.Exp(-6.74 * F_0) * (Math.PI * H * L / (4.0 * F_f) - 0.02) - A * (1.0 - Math.Exp(-4.5 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / (R_x + R_xsh)));

            fi[3] = H / L - Math.Exp(1.0 - 2.0 * R_yd / (Q + 2.0 * G_g + G_z));
            return;
            //scaling
            fi[0] /= R_x;
            fi[1] /= H * L * Math.PI / (4 * F_f);
            fi[2] /= A * (1.0 - Math.Exp(-4.5 * (2.0 * R_yd + Q + 2.0 * G_g + G_z) / (R_x + R_xsh)));
            fi[3] /= H / L;
        }
       
        static double epsg = 1e-5;
        static double epsf = 0;
        static double epsx = 0;
        static int maxits = 0;
        //
        static void minlm(ref double[] x)
        {

           
            alglib.minlmstate state;
            alglib.minlmreport rep;
            alglib.minlmcreatev(4, x, 0.0001, out state);
            alglib.minlmsetcond(state, epsg, epsf, epsx, maxits);
            alglib.minlmsetscale(state, new double[] { 10, 10, 10, 10000 });
            alglib.minlmoptimize(state, f, null, null);
            alglib.minlmresults(state, out x, out rep);
            System.Console.WriteLine("{0}", rep.terminationtype); // EXPECTED: 4
        }
        static void Main()
        {




            double[] x = { 50, 50, 6, 10000 };

            double[] fi = { 0, 0, 0, 0 };

            minlm(ref x);
            f(x, fi, null);
           
            System.Console.WriteLine(" X : {0}, {1}, {2}, {3}", x[0], x[1], x[2], x[3]);
            System.Console.WriteLine(" f(X) : {0}, {1}, {2}, {3}", fi[0], fi[1], fi[2], fi[3]);
            System.Console.ReadLine();
        }
    }
}

I would appreciate any suggestions


Top
 Profile  
 
 Post subject: Re: nonlinear equations system
PostPosted: Wed Jun 08, 2011 8:14 pm 
Offline
Site Admin

Joined: Fri May 07, 2010 7:06 am
Posts: 900
Quote:
hey! f(x)[0] is not 0! but the solution is almost exact!

Yes, I just noticed it. But it leads us to the question - are you sure that you've correctly converted your equation from Mathcad to C#? f0(x)=42 is too large to be considered "approximately zero" by Mathcad.

It can be the case that function your code actually calculates has minimum at different location than original one. And, unlike original problem, your incorrectly typed problem has no solution.

Can you focus on that possibility? and I will try to debug solver.


Top
 Profile  
 
Display posts from previous:  Sort by  
Post new topic Reply to topic  [ 11 posts ]  Go to page 1, 2  Next

All times are UTC


Who is online

Users browsing this forum: No registered users and 4 guests


You cannot post new topics in this forum
You cannot reply to topics in this forum
You cannot edit your posts in this forum
You cannot delete your posts in this forum
You cannot post attachments in this forum

Search for:
Jump to:  
cron
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group