Problema matemático complejo

[DarkDudae]

Buenas a tod@s:

¿Qué tal lleváis los sistemas de ecuaciones no lineales? Espero que mejor que yo. A ver si me podéis echar una mano:


Dado un vector con N valores, que podemos llamar VectorValores. Por otro lado, Tengo una función que genera a su vez otro vector de otros N valores. Llamémoslo VectorAjuste.

La función que generará los velores del VectorAjuste es la siguiente:

Código Delphi [-]
  Exp1:=(-1)*((Power(i-xk,2)/(2*Power(sigmaK,2))));
  Exp2:=(-1)*((Power(i-xLa,2)/(2*Power(sigmaLa,2))));
  Result:=A*EXP(Exp1)+B*EXP(Exp2);

En la función anterior, "i" será el índice del vector de N valores. El resto de parámetros: A, B, xk, sigmaK, xLa y sigmaLa serán incógnitas.

Ahora viene la parte "divertida":

Tengo que encontrar las 6 incógnitas citadas ( A, B, xk, sigmaK, xLa y sigmaLa) de forma que el Sumatorio de i hasta N de (VectorValores[i]-VectorAjuste[i]) al cuadrado sea mínimo:

Código Delphi [-]
    Resultado:=0;
    for i:=0 to N-1 do
      begin
        Resultado:=Resultado+Power(VectorValores[i]-VectorAjuste[i],2));
      end;

Es decir, el valor "Resultado", una vez calculadas las 6 incógnitas, tiene que ser mínimo.

En realidad, esto es un problema de "root-finding", que por ejemplo Excel puede resolver mediante la herramienta Solver que tiene incorporada. Sé que esta herramienta usa el algoritmo de Newton Raphson, pero sinceramente no sé cómo aplicarlo en un sistema tan complejo.

A la herramienta solver de Excel hay que proporcionarles unos valores "aproximados" para que empiece a buscar los valores para que cumplan el criterio. En el caso correspondiente, los valores son los siguientes:

A xK σK B xLa σLa
250 750 20 200 775 20

Sé cómo implementar un sistema de resolución de ecuaciones "Simple" mediante Newton Raphson, como por ejemplo resolver las siguientes ecuaciones:
f1(x,y) = x^2 + y^2 - 4 = 0
f2(x,y) = x^2 - y - 1 = 0

Pero claro, esto dista mucho del problema resuelto. Lo puedo hacer por fuerza bruta dejando como constantes los valores de sigma, pero es bastante lento de calcular. Sé que siempre tengo la opción de usar la macro de excel para invocar a la herramienta "solver", pero no me apetece tener que depender de herramientas de terceros.

Cualquier ayuda y/o pista será más que bien recibida.

[fenixariel]

Un buen libro de Metodos Numericos?

Yo tengo uno : 'Metodos Numericos con Matlab', pasarlo a pascal muy facil. Me parece muy muy bueno el libro, me ha servido mucho pero mucho.


Si es matematica, es matematica.

Creo que primero solucionar las ecuaciones en papel, claro todo literal.


Saludos.

[[nlsgarcia]]

DarkDudae,

Empezado por DarkDudae

¿Qué tal lleváis los sistemas de ecuaciones no lineales?

Revisa estos links:

Métodos numéricos de resolución de ecuaciones no lineales
http://ocw.upm.es/matematica-aplicad...s/EcsNoLin.pdf

Introdución a los Métodos Númericos
http://www.tec-digital.itcr.ac.cr/re...sNumericos.pdf

Numerical Recipes Books On-Line
http://www.nr.com/oldverswitcher.html

Métodos Numéricos en Ingeniería
http://disi.unal.edu.co/~lctorress/MetNum/LiMetNu2.pdf

Adicionalmente te sugiero revisar este libro: Métodos Numéricos para Ingenieros de Steven C. Chapra y Raymond P. Canale, Editorial McGraw-Hill.

Espero sea útil

Nelson.

[DarkDudae]

Gracias por los enlaces.

[DarkDudae]

Escribo nuevamente para comentaros que ya he podido solucionar el problema. Al final lo he conseguido mediante el método de Levenberg-Maxquadrat. En las librerías de ALGLIB podemos encontrar una implementación del mismo.

La parte más "compleja" es que hay que calcular las derivadas de cada una de las incógnitas para mejorar el proceso de convergencia. Pero bueno, lo dejo aquí por si le puede ser de ayuda a alguien en un futuro por si tiene que minimizar funciones complejas.

 Memo1.Clear;
    //Tenemos 4 incógnitas: A, Xk, B, XLa
    N := 4;

    SetLength(S, N);

    //S: Nuestra aproximación inicial que sabemos no excesivamente alejada de la realidad
    S[0]:=250; S[1]:=750; S[2]:=200; S[3]:=775;

    M := 120;
    SetLength(X1, M);
    SetLength(Y1, M);
    I:=0;
    while I<=M-1 do
    begin
        X1[i] := 700+I;
        //Leemos los 120 valores de un memo (mValores) o cualquier stream
        Y1[i] := strtofloat(mValores.Lines[i]);
        Inc(I);
    end;

    //Ahora S almacena los puntos de inicio, mientras que X e Y almacenan los puntos que serán "suavizados y/o minimizados"
    MinLMCreateFJ(N, M, S, State);
    MinLMSetCond(State, 0.0, 0.0, 0.001, 0);
    while MinLMIteration(State) do
    begin

        if State.NeedF then
        begin
            State.F := 0;
        end;
        I:=0;
        while I<=M-1 do
        begin

            //
            // "A" guardada en State.X[0]
            // "Xk" - State.X[1]
            // "B" - State.X[2]
            // "XLa" - State.X[3]
            //

            FI := Y1[i]-State.X[0]*exp((-1)*((AP_Sqr(X1[i]-State.X[1])/(2*AP_Sqr(6.3)))))-State.X[2]*exp((-1)*((AP_Sqr(X1[i]-State.X[3])/(2*AP_Sqr(6.3)))));

            if State.NeedF then
            begin
               //Sumatorio aquí
                State.F := State.F+AP_Sqr(FI);
            end;
            if State.NeedFiJ then
            begin
                //
                // Fi
                //
                State.Fi[i] := FI;

                //
                // dFi/dA = -exp^(-((x-Xk)^2)/2s^2)
                //
                State.J[I,0] := (-1)*exp((-1)*((AP_Sqr(X1[i]-State.X[1])/(2*AP_Sqr(6.3)))));
                //
                // dFi/dXk = (A*(xK-x)*exp(-(((x-Xk)^2)/(2*s^2))))))/(6.3^2))

                //
                State.J[I,1] := (State.X[0]*(State.X[1]-X1[i])*exp((-1)*((AP_Sqr(X1[i]-State.X[1])/(2*AP_Sqr(6.3))))))/(AP_Sqr(6.3));

                //
                // dFi/dB = -exp^(-((x-XL)^2)/2s^2) 
                //
                State.J[I,2] := (-1)*exp((-1)*((AP_Sqr(X1[i]-State.X[3])/(2*AP_Sqr(6.3)))));

                //
                // dFi/dXLa = (B*(xL-x)*exp(-(((x-xL)^2)/(2*s^2))))))/(6.3^2))
                //
                State.J[I,3] := (State.X[2]*(State.X[3]-X1[i])*exp((-1)*((AP_Sqr(X1[i]-State.X[3])/(2*AP_Sqr(6.3))))))/(AP_Sqr(6.3));


            end;
            Inc(I);
        end;
    end;
    MinLMResults(State, S, Rep);

    //
    // output results
    //
    Memo1.Lines.Add(Format('A = %4.2f'#13#10'',[
        S[0]]));
    Memo1.Lines.Add(Format('Xk = %4.2f'#13#10'',[
        S[1]]));
    Memo1.Lines.Add(Format('B = %4.2f'#13#10'',[
        S[2]]));
    Memo1.Lines.Add(Format('XLa = %4.2f'#13#10'',[
        S[3]]));
    Memo1.Lines.Add(Format('Tipo de terminación = %0d (debería ser 2 - parando cuando sea suficientemente pequeño)'#13#10'',[
        Rep.TerminationType]));