Examples of LUdcmp

com.nr.la.LUdcmp
LU decomposition - PTC Iterative Improvement of a Solution to Linear Equations Copyright (C) Numerical Recipes Software 1986-2007 Java translation Copyright (C) Huang Wen Hui 2012 @author hwh

Examples of com.nr.la.LUdcmp

        sum += (rbf[i][j] = fn.rbf(rad(pts[i],pts[j])));
      }
      if (norm) rhs[i] = sum*vals[i];
      else rhs[i] = vals[i];
    }
    LUdcmp lu = new LUdcmp(rbf);   // Solve the set of linear equations
    lu.solve(rhs,w);
  }

View Full Code Here

Examples of com.nr.la.LUdcmp

      }
      v[i][npt] = v[npt][i] = 1.;
    }
    v[npt][npt] = y[npt] = 0.;
    if (err!=null) for (i=0;i<npt;i++) v[i][i] -= SQR(err[i]);
    vi = new LUdcmp(v);
    vi.solve(y,yvi);
  }

View Full Code Here

Examples of com.nr.la.LUdcmp

    double[] denom = new double[n+1];
    for (j=0;j<n;j++) {
      y[j]=cof[n+j+1];
      for (k=0;k<n;k++) q[j][k]=cof[j-k+n];
    }
    LUdcmp lu = new LUdcmp(q);
    lu.solve(y,x);
    for (j=0;j<4;j++) lu.mprove(y,x);
    for (k=0;k<n;k++) {
      for (sum=cof[k+1],j=0;j<=k;j++) sum -= x[j]*cof[k-j];
      y[k]=sum;
    }
    num[0] = cof[0];

View Full Code Here

Examples of com.nr.la.LUdcmp

    double h=htot/nstep;
    for (int i=0;i<n;i++) {
      for (int j=0;j<n;j++) a[i][j] = -dfdy[i][j];
      a[i][i] += 1.0/h;
    }
    LUdcmp alu = new LUdcmp(a);
    double xnew=x+h;
    derivs.derivs(xnew,y,del);
    for (int i=0;i<n;i++)
      ytemp[i]=y[i];
    alu.solve(del,del);
    if (dense && nstep==k+1) {
      ipt.val++;
      for (int i=0;i<n;i++)
        fsave[ipt.val][i]=del[i];
    }
    for (int nn=1;nn<nstep;nn++) {
      for (int i=0;i<n;i++)
        ytemp[i] += del[i];
      xnew += h;
      derivs.derivs(xnew,ytemp,yend);
      if (nn ==1 && k<=1) {
        double del1=0.0;
        for (int i=0;i<n;i++)
          del1 += SQR(del[i]/scale[i]);
        del1=sqrt(del1);
        derivs.derivs(x+h,ytemp,dytemp);
        for (int i=0;i<n;i++)
          del[i]=dytemp[i]-del[i]/h;
        alu.solve(del,del);
        double del2=0.0;
        for (int i=0;i<n;i++)
          del2 += SQR(del[i]/scale[i]);
        del2=sqrt(del2);
        theta=del2/max(1.0,del1);
        if (theta > 1.0)
          return false;
      }
      alu.solve(yend,del);
      if (dense && nn >= nstep-k-1) {
        ipt.val++;
        for (int i=0;i<n;i++)
          fsave[ipt.val][i]=del[i];
      }

View Full Code Here

Examples of com.nr.la.LUdcmp

    int i;
    for (i=0;i<n;i++) {
      for (int j=0;j<n;j++) a[i][j] = -dfdy[i][j];
      a[i][i] += 1.0/(gam*h);
    }
    LUdcmp alu = new LUdcmp(a);
    for (i=0;i<n;i++)
      ytemp[i]=dydx[i]+h*d1*dfdx[i];
    alu.solve(ytemp,k1);
    for (i=0;i<n;i++)
      ytemp[i]=y[i]+a21*k1[i];
    derivs.derivs(x+c2*h,ytemp,dydxnew);
    for (i=0;i<n;i++)
      ytemp[i]=dydxnew[i]+h*d2*dfdx[i]+c21*k1[i]/h;
    alu.solve(ytemp,k2);
    for (i=0;i<n;i++)
      ytemp[i]=y[i]+a31*k1[i]+a32*k2[i];
    derivs.derivs(x+c3*h,ytemp,dydxnew);
    for (i=0;i<n;i++)
      ytemp[i]=dydxnew[i]+h*d3*dfdx[i]+(c31*k1[i]+c32*k2[i])/h;
    alu.solve(ytemp,k3);
    for (i=0;i<n;i++)
      ytemp[i]=y[i]+a41*k1[i]+a42*k2[i]+a43*k3[i];
    derivs.derivs(x+c4*h,ytemp,dydxnew);
    for (i=0;i<n;i++)
      ytemp[i]=dydxnew[i]+h*d4*dfdx[i]+(c41*k1[i]+c42*k2[i]+c43*k3[i])/h;
    alu.solve(ytemp,k4);
    for (i=0;i<n;i++)
      ytemp[i]=y[i]+a51*k1[i]+a52*k2[i]+a53*k3[i]+a54*k4[i];
    double xph=x+h;
    derivs.derivs(xph,ytemp,dydxnew);
    for (i=0;i<n;i++)
      k6[i]=dydxnew[i]+(c51*k1[i]+c52*k2[i]+c53*k3[i]+c54*k4[i])/h;
    alu.solve(k6,k5);
    for (i=0;i<n;i++)
      ytemp[i] += k5[i];
    derivs.derivs(xph,ytemp,dydxnew);
    for (i=0;i<n;i++)
      k6[i]=dydxnew[i]+(c61*k1[i]+c62*k2[i]+c63*k3[i]+c64*k4[i]+c65*k5[i])/h;
    alu.solve(k6,yerr);
    for (i=0;i<n;i++)
      yout[i]=ytemp[i]+yerr[i];
  }

View Full Code Here

Examples of com.nr.la.LUdcmp

    final int N=40;
    
    double[] g = new double[N];
    double[][] a = new double[N][N];
    new Quad_matrix(a);
    LUdcmp alu = new LUdcmp(a);
    for (int j=0;j<N;j++)
      g[j]=sin(j*PI/(N-1));
    alu.solve(g,g);
    for (int j=0;j<N;j++) {
      double x=j*PI/(N-1);
      System.out.printf("%d  %f  %f\n", (j+1), x, g[j]);
    }
  }

View Full Code Here

Examples of com.nr.la.LUdcmp

      for (int j=0;j<n;j++) {
        omk[i][j]=((i == j) ? 1 : 0)-ak(t[i],t[j])*w[j];
      }
      f[i]=g(t[i]);
    }
    LUdcmp alu = new LUdcmp(omk);
    alu.solve(f,f);
  }

View Full Code Here

Examples of com.nr.la.LUdcmp

          else
            a[k][l] = -0.5 * h * ak(k, l, t[i], t[i]);
        }
        b[k] = sum;
      }
      LUdcmp alu = new LUdcmp(a);
      alu.solve(b, b);
      for (int k = 0; k < m; k++)
        f[k][i] = b[k];
    }
  }

View Full Code Here

Examples of com.nr.la.LUdcmp

      for (k=1;k<=nr;k++) sum += pow(k,ipj);
      for (k=1;k<=nl;k++) sum += pow(-k,ipj);
      mm=min(ipj,2*m-ipj);
      for (imj = -mm;imj<=mm;imj+=2) a[(ipj+imj)/2][(ipj-imj)/2]=sum;
    }
    LUdcmp alud = new LUdcmp(a);
    for (j=0;j<m+1;j++) b[j]=0.0;
    b[ld]=1.0;
    alud.solve(b,b);
    for (kk=0;kk<np;kk++) c[kk]=0.0;
    for (k = -nl;k<=nr;k++) {
      sum=b[0];
      fac=1.0;
      for (mm=1;mm<=m;mm++) sum += b[mm]*(fac *= k);

View Full Code Here

Examples of com.nr.la.LUdcmp


    


    // Test ludcmp
    System.out.println("Testing ludcmp");
    LUdcmp alu = new LUdcmp(a);
    double[][] x = new double[b.length][b[0].length];
    alu.solve(b,x);
    sbeps=5.e-15;
    double maxel = maxel(matsub(matmul(a,x),b));
    localflag = maxel > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** ludcmp: test of solution vector failed: "+ maxel);
      
    }


    // Test inverse
    double[][] ainv; // = buildMatrix(a); not need it
    ainv = alu.inverse();
    sbeps = 5.e-15;
    localflag = maxel(matsub(matmul(ainv,a),ident(a.length,1.))) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** ludcmp: Test of inverse failed");
      
    }


    // Test determinant
    LUdcmp ainvlu = new LUdcmp(ainv);
    sbeps = 5.e-15;
    localflag = (alu.det()*ainvlu.det()-1.) > sbeps;
    globalflag = globalflag || localflag;
    if (localflag) {
      fail("*** ludcmp: Test of determinant failed");
      
    }

View Full Code Here

0 1

TOP

All source code are property of their respective owners. Java is a trademark of Sun Microsystems, Inc and owned by ORACLE Inc. Contact coftware#gmail.com.