Package com.nr.ode

Source Code of com.nr.ode.StepperRoss

package com.nr.ode;

import static com.nr.NRUtil.*;
import static java.lang.Math.*;
import org.netlib.util.*;

import com.nr.la.LUdcmp;

/**
* fourth-order stiffly stable Rosenbrock stepper
* Copyright (C) Numerical Recipes Software 1986-2007
* Java translation Copyright (C) Huang Wen Hui 2012
*
* @author hwh
*
*/
public class StepperRoss extends StepperBase{
  final static  double c2=0.386;
  final static  double c3=0.21;
  final static  double c4=0.63;
  final static  double bet2p=0.0317;
  final static  double bet3p=0.0635;
  final static  double bet4p=0.3438;
  final static  double d1= 0.2500000000000000e+00;
  final static  double d2=-0.1043000000000000e+00;
  final static  double d3= 0.1035000000000000e+00;
  final static  double d4=-0.3620000000000023e-01;
  final static  double a21= 0.1544000000000000e+01;
  final static  double a31= 0.9466785280815826e+00;
  final static  double a32= 0.2557011698983284e+00;
  final static  double a41= 0.3314825187068521e+01;
  final static  double a42= 0.2896124015972201e+01;
  final static  double a43= 0.9986419139977817e+00;
  final static  double a51= 0.1221224509226641e+01;
  final static  double a52= 0.6019134481288629e+01;
  final static  double a53= 0.1253708332932087e+02;
  final static  double a54=-0.6878860361058950e+00;
  final static  double c21=-0.5668800000000000e+01;
  final static  double c31=-0.2430093356833875e+01;
  final static  double c32=-0.2063599157091915e+00;
  final static  double c41=-0.1073529058151375e+00;
  final static  double c42=-0.9594562251023355e+01;
  final static  double c43=-0.2047028614809616e+02;
  final static  double c51= 0.7496443313967647e+01;
  final static  double c52=-0.1024680431464352e+02;
  final static  double c53=-0.3399990352819905e+02;
  final static  double c54= 0.1170890893206160e+02;
  final static  double c61= 0.8083246795921522e+01;
  final static  double c62=-0.7981132988064893e+01;
  final static  double c63=-0.3152159432874371e+02;
  final static  double c64= 0.1631930543123136e+02;
  final static  double c65=-0.6058818238834054e+01;
  final static  double gam= 0.2500000000000000e+00;
  final static  double d21= 0.1012623508344586e+02;
  final static  double d22=-0.7487995877610167e+01;
  final static  double d23=-0.3480091861555747e+02;
  final static  double d24=-0.7992771707568823e+01;
  final static  double d25= 0.1025137723295662e+01;
  final static  double d31=-0.6762803392801253e+00;
  final static  double d32= 0.6087714651680015e+01;
  final static  double d33= 0.1643084320892478e+02;
  final static  double d34= 0.2476722511418386e+02;
  final static  double d35=-0.6594389125716872e+01;
  double[][] dfdy;
  double[] dfdx;
  double[] k1,k2,k3,k4,k5,k6;
  double[] cont1,cont2,cont3,cont4;
  double[][] a;
  Controller con = new Controller();
 
  public static class Controller {
    double hnext;
    boolean reject;
    boolean first_step;
    double errold;
    double hold;
   
    public Controller(){
      reject=false;
      first_step=true;
    }
   
    public boolean success(final double err, final doubleW h) {
      final double safe=0.9,fac1=5.0,fac2=1.0/6.0;
      double fac=max(fac2,min(fac1,pow(err,0.25)/safe));
      double hnew=h.val/fac;
      if (err <= 1.0) {
        if (!first_step) {
          double facpred=(hold/h.val)*pow(err*err/errold,0.25)/safe;
          facpred=max(fac2,min(fac1,facpred));
          fac=max(fac,facpred);
          hnew=h.val/fac;
        }
        first_step=false;
        hold=h.val;
        errold=max(0.01,err);
        if (reject)
          hnew=(h.val >= 0.0 ? min(hnew,h.val) : max(hnew,h.val));
        hnext=hnew;
        reject=false;
        return true;
      } else {
        h.val=hnew;
        reject=true;
        return false;
      }
    }
  }
 
  public StepperRoss(){
   
  }
 
  public StepperRoss(final double[] yy, final double[] dydxx, final double xx,
      final double atoll,final double rtoll, final boolean dense){
    setParam(yy,dydxx,xx,atoll,rtoll,dense);
  }
 
  public void setParam(final double[] yy, final double[] dydxx, final double xx,
    final double atoll,final double rtoll, final boolean dense) {
    super.setParam(yy,dydxx,xx,atoll,rtoll,dense);
    dfdy = new double[n][n];
    dfdx = new double[n];
    k1 = new double[n];
    k2 = new double[n];
    k3 = new double[n];
    k4 = new double[n];
    k5 = new double[n];
    k6 = new double[n];
    cont1 = new double[n];
    cont2 = new double[n];
    cont3 = new double[n];
    cont4 = new double[n];
    a= new double[n][n];
    EPS=DBL_EPSILON;
  }
 
  public void step(final double htry,final DerivativeInf derivs) {
    double[] dydxnew = new double[n];
    doubleW h= new doubleW(htry);
    derivs.jacobian(x,y,dfdx,dfdy);
    for (;;) {
      dy(h.val,derivs);
      double err=error();
      if (con.success(err,h)) break;
      if (abs(h.val) <= abs(x)*EPS)
        throw new IllegalArgumentException("stepsize underflow in StepperRoss");
    }
    derivs.derivs(x+h.val,yout,dydxnew);
    if (dense)
      prepare_dense(h.val,dydxnew);
    System.arraycopy(dydxnew, 0, dydx, 0, dydxnew.length);
    System.arraycopy(yout, 0, y, 0, y.length);
    //dydx=dydxnew;
    //y=yout;
    xold=x;
    x += (hdid=h.val);
    hnext=con.hnext;
  }
 
  public void dy(final double h,final DerivativeInf derivs) {
    double[] ytemp = new double[n],dydxnew = new double[n];
    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];
  }
 
  public void prepare_dense(final double h,final double[] dydxnew) {
    for (int i=0;i<n;i++) {
      cont1[i]=y[i];
      cont2[i]=yout[i];
      cont3[i]=d21*k1[i]+d22*k2[i]+d23*k3[i]+d24*k4[i]+d25*k5[i];
      cont4[i]=d31*k1[i]+d32*k2[i]+d33*k3[i]+d34*k4[i]+d35*k5[i];
    }
  }
 
  public double dense_out(final int i,final double x,final double h) {
    double s=(x-xold)/h;
    double s1=1.0-s;
    return cont1[i]*s1+s*(cont2[i]+s1*(cont3[i]+s*cont4[i]));
  }
 
  public double error() {
    double err=0.0,sk;
    for (int i=0;i<n;i++) {
      sk=atol+rtol*max(abs(y[i]),abs(yout[i]));
      err += SQR(yerr[i]/sk);
    }
    return sqrt(err/n);
  }
}
TOP

Related Classes of com.nr.ode.StepperRoss

TOP
Copyright © 2018 www.massapi.com. All rights reserved.
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.