001    package numal.highprecision.test;
002    
003    import java.text.DecimalFormat;
004    import numal.highprecision.*;
005    
006    /** Determine the parameters p[1], p[2], p[3] of the best fit of the function
007     * g(x[i],p)=p[1]+p[2]*exp(p[3]*x[i]) to data readings 
008     * y[]={127.0, 151.0, 379.0, 421.0, 460.0, 426.0} where
009     * x[]={-5, -3, -1, 1, 3, 5}.
010     * The components of the residual vector are  
011     *    p[1]+p[2]*exp(p[3]*x[i])-y[i]. 
012     * The elements of the Jacobian matrix are 
013     *    dg[i]/dp[1]=1, dg[i]/dp[2]=exp(p[3]*x[i]), dg[i]/dp[3]=p2[2]*x[i]*exp(p[3]*x[i]).
014     *
015     * For a successful test the main() method must output the following values :
016     * {@code
017        Parameters:
018          5.23225E2   -1.56839E2   -1.99778E-1
019        
020        OUT:
021            7.22078E7
022            1.15716E2
023            1.72802E-3
024            1.65459E2
025            2.30000E1
026            2.20000E1
027            0.00000E0
028        
029        Last residual vector:
030          -29.6  86.6  -47.3  -26.2  -22.9  39.5
031        }
032     *
033     */
034    public class Test_marquardt extends Object
035                                implements AP_marquardt_methods {
036            
037      static double x[] = new double[7];
038      static double y[] = new double[7];
039    
040            public static void main(String args[]) {
041    
042        double in[] = new double[7];
043        double out[] = new double[8];
044        double rv[] = new double[7];
045        double par[] = new double[4];
046        double jjinv[][] =  new double[4][4];
047        
048        Test_marquardt testmarquardt = new Test_marquardt();
049        DecimalFormat fiveDigit = new DecimalFormat("0.00000E0");
050        DecimalFormat oneDigit = new DecimalFormat("0.0");
051        in[0]=1.0e-6;  in[3]=1.0e-4;  in[4]=1.0e-1;  in[5]=75.0;
052        in[6]=1.0e-2;
053        x[1] = -5.0;  x[2] = -3.0;  x[3] = -1.0;  x[4]=1.0;
054        x[5]=3.0;  x[6]=5.0;
055        y[1]=127.0;  y[2]=151.0;  y[3]=379.0;  y[4]=421.0;
056        y[5]=460.0;  y[6]=426.0;
057        par[1]=580.0;  par[2] = -180.0;  par[3] = -0.160;
058        Analytic_problems.marquardt(6,3,par,rv,jjinv,
059                                    testmarquardt,in,out);
060        System.out.println("Parameters:\n  " +
061          fiveDigit.format(par[1]) + "   " +
062          fiveDigit.format(par[2]) + "   " +
063          fiveDigit.format(par[3]) + "\n\nOUT:\n    " +
064          fiveDigit.format(out[7]) + "\n    " +
065          fiveDigit.format(out[2]) + "\n    " +
066          fiveDigit.format(out[6]) + "\n    " +
067          fiveDigit.format(out[3]) + "\n    " +
068          fiveDigit.format(out[4]) + "\n    " +
069          fiveDigit.format(out[5]) + "\n    " +
070          fiveDigit.format(out[1]) +
071          "\n\nLast residual vector:\n  " +
072          oneDigit.format(rv[1]) + "  " + oneDigit.format(rv[2]) +
073          "  " + oneDigit.format(rv[3]) + "  " + 
074          oneDigit.format(rv[4]) + "  " + oneDigit.format(rv[5]) +
075          "  " + oneDigit.format(rv[6]));
076            }
077    
078    
079      public boolean funct(int m, int n,double par[], double rv[])
080      {
081        int i;
082    
083        for (i=1; i<=m; i++) {
084          if (par[3]*x[i] > 680.0) return false;
085          rv[i]=par[1]+par[2]*Math.exp(par[3]*x[i])-y[i];
086        }
087        return true;
088      }
089      
090      
091      public void jacobian(int m, int n, double par[],
092                           double rv[], double jac[][])
093      {
094        int i;
095        double ex;
096    
097        for (i=1; i<=m; i++) {
098          jac[i][1]=1.0;
099          jac[i][2]=ex=Math.exp(par[3]*x[i]);
100          jac[i][3]=x[i]*par[2]*ex;
101        }
102      }
103    }