ptolemy.domains.algebraic.kernel

Class AlgebraicLoopDirector.Homotopy

java.lang.Object

ptolemy.domains.algebraic.kernel.AlgebraicLoopDirector.AlgebraicLoopSolver

ptolemy.domains.algebraic.kernel.AlgebraicLoopDirector.Homotopy

Enclosing class:

AlgebraicLoopDirector

class AlgebraicLoopDirector.Homotopy
extends AlgebraicLoopDirector.AlgebraicLoopSolver

Class for solving algebraic loops using a homotopy method.

Author:

Michael Wetter

Field Summary

Fields

Modifier and Type	Field and Description
protected double	_acfac Acceleration factor for step length control
protected double	_angmax Maximal angle
protected double[][]	_b Matrix b used in Newton algorithm.
protected double	_c1 Value c1 used in Newton algorithm.
protected double	_c2 Value c2 used in Newton algorithm.
protected double	_cdmax Maximum for condition estimate
protected double	_ctmax Maximum contraction rate in corrector step, 0 < ctmax < 1.
protected double[]	_deltaX Step size for finite difference approximation
protected double	_dmax Maximal norm for H
protected double	_dmin Minimal norm for H
protected double	_h Initial step size
protected double	_hmax Maximal step size
protected double	_hmin Minimal step size
protected double	_hmn Minimal Newton step size
protected int	_n1 Number of independent variables including the homotopy factor
protected double[][]	_q Matrix q used in Newton algorithm.
protected double[]	_r Result of Newton step
protected double[]	_t Tangent vector to homotopy curve.
protected boolean	_test Test for step length in Newton algorithm
protected double[]	_xIni Initial guess
protected double[]	_y Current guess of solution

Fields inherited from class ptolemy.domains.algebraic.kernel.AlgebraicLoopDirector.AlgebraicLoopSolver

_converged, _iterationCount, _maxIterations, _tolerance, _variableNames

Constructor Summary

Constructors

Constructor and Description
Homotopy(java.lang.String[] variableNames, double[] tolerance, int maxIterations) Construct an algebraic loop solver.

Method Summary

Modifier and Type	Method and Description
(package private) double	_decomp() Conduct a QR decomposition.
protected double	_getOrientation() Return the direction in which the curve will be traversed.
(package private) void	_givens(int l1, int l2, int l3) Perform a Givens rotation.
(package private) double[][]	_jac(double[] x, double h) Evaluate the transpose of the Jacobian of H(x, lambda) = x - lambda F(x) by using forward differences.
protected double	_l2norm(double[] x) Return the L2 norm.
(package private) double[]	_map(double[] x) Compute y = H(x), where H(x)=0 is curve to be traced.
protected void	_newton(double[] u, double[] v, double[] w) Conduct a Newton step.
protected void	_updateQB(double[] w, double angmax) Update _q and _b arrays.
void	solve(double[] xIni) This method solves u - lambda F(u) = 0 with initial values u=0 and lambda=0.

Methods inherited from class ptolemy.domains.algebraic.kernel.AlgebraicLoopDirector.AlgebraicLoopSolver

_didConverge, converged, getIterationCount

Methods inherited from class java.lang.Object

clone, equals, finalize, getClass, hashCode, notify, notifyAll, toString, wait, wait, wait

Field Detail

_ctmax

protected double _ctmax

Maximum contraction rate in corrector step, 0 < ctmax < 1.

_dmax

protected double _dmax

Maximal norm for H

_dmin

protected double _dmin

Minimal norm for H

_hmax

protected double _hmax

Maximal step size

_hmin

protected double _hmin

Minimal step size

_hmn

protected double _hmn

Minimal Newton step size

_h

protected double _h

Initial step size

_n1

protected int _n1

Number of independent variables including the homotopy factor

_cdmax

protected double _cdmax

Maximum for condition estimate

_angmax

protected double _angmax

Maximal angle

_acfac

protected double _acfac

Acceleration factor for step length control

_b

protected double[][] _b

Matrix b used in Newton algorithm.

_q

protected double[][] _q

Matrix q used in Newton algorithm.

_t

protected double[] _t

Tangent vector to homotopy curve.

_r

protected double[] _r

Result of Newton step

_xIni

protected double[] _xIni

Initial guess

_y

protected double[] _y

Current guess of solution

_test

protected boolean _test

Test for step length in Newton algorithm

_c1

protected double _c1

Value c1 used in Newton algorithm.

_c2

protected double _c2

Value c2 used in Newton algorithm.

_deltaX

protected double[] _deltaX

Step size for finite difference approximation

Constructor Detail

Homotopy

public Homotopy(java.lang.String[] variableNames,
                double[] tolerance,
                int maxIterations)
         throws IllegalActionException

Construct an algebraic loop solver.

Parameters:

variableNames - Names of each break variable.

tolerance - Tolerance for each variable.

maxIterations - Maximum number of iterations.

Throws:

IllegalActionException

Method Detail

solve

public void solve(double[] xIni)
           throws IllegalActionException

This method solves u - lambda F(u) = 0 with initial values u=0 and lambda=0. The solution, for lambda=1, is a fixed point of F : Re^n -> Re^n. To allow a non-zero initial guess for the loop function, the problem is reformulated as H(s, lambda, x0) = s - lambda (g(s+x0)-x0), where lambda has an initial value of 0 and is successively increased to 1, s is a coordinate transformation defined so that x = s+x0, where x0 is the initial iterate. The Jacobian is approximated at the start of the solution using finite differences, and then updated using Broyden's method. The method starts with computing the Jacobian H'(x)=A, where x=(s, lambda). Next, it computes the tangent vector t = t(A). It then conducts a predictor Euler step u = x+h*t. After computing a perturbation vector pv, it corrects the iterates using v = u - A^+ (H(u)-pv), where A^+ is the Moore-Penrose inverse of A.

Specified by:

solve in class AlgebraicLoopDirector.AlgebraicLoopSolver

Parameters:

xIni - Array with the initial values of the variables, which will be replaced by this method with the obtained solution.

Throws:

IllegalActionException - If the prefire() method returns false having previously returned true in the same iteration, or if the prefire() or fire() method of the actor throws it, or if evaluating the function yields a value that is not a double, or if the solver fails to find a solution.

_jac

double[][] _jac(double[] x,
                double h)
         throws IllegalActionException

Evaluate the transpose of the Jacobian of H(x, lambda) = x - lambda F(x) by using forward differences.

Parameters:

x - Point at which the Jacobian is approximated.

h - Step size.

Returns:

The transpose of the Jacobian.

Throws:

IllegalActionException - If the prefire() method returns false having previously returned true in the same iteration, or if the prefire() or fire() method of the actor throws it, or if evaluating the function yields a value that is not a double.

_map

double[] _map(double[] x)
       throws IllegalActionException

Compute y = H(x), where H(x)=0 is curve to be traced.

Parameters:

x - Independent variable where the last element is the homotopy factor.

Returns:

y = H(x).

Throws:

IllegalActionException - If the prefire() method returns false having previously returned true in the same iteration, or if the prefire() or fire() method of the actor throws it, or if evaluating the function yields a value that is not a double.

_getOrientation

protected double _getOrientation()

Return the direction in which the curve will be traversed.

Returns:

Direction in which the curve will be traversed.

_givens

void _givens(int l1,
             int l2,
             int l3)

Perform a Givens rotation. This method performs a Givens rotation on _b and _q. Prior to calling this method, _c1 and _c2 need to be set. The method then uses _c1 and _c2, and sets them to new values. A method that calls _givens then need to use _c1 and _c2. This was needed as Java passes double by value and not be reference.

Parameters:

l1 - Coordinate to be acted upon.

l2 - Coordinate to be acted upon.

l3 - Coordinate to be acted upon.

_decomp

double _decomp()

Conduct a QR decomposition. A QR decomposition for _b is stored in _q and _b by using Givens rotation on _b an _q until _b is upper triangular. A very coarse condition estimate is returned.

Returns:

A very coarse condition estimate.

_newton

protected void _newton(double[] u,
                       double[] v,
                       double[] w)
                throws IllegalActionException

Conduct a Newton step. Conduct a Newton step v = u - A^+ where A is approximated by H'. The matrix A^+ is the Moore-Penrose inverse of A. This method uses perturbations to stabilize the method and performs tests on the residuals and the contractions.

Parameters:

u -

w - This argument is changed by this function.

Throws:

IllegalActionException - If the prefire() method returns false having previously returned true in the same iteration, or if the prefire() or fire() method of the actor throws it, or if evaluating the function yields a value that is not a double.

_updateQB

protected void _updateQB(double[] w,
                         double angmax)

Update _q and _b arrays. This method updates the _q and _b arrays using QR decomposition.

_l2norm

protected double _l2norm(double[] x)

Return the L2 norm.

Parameters:

x - Argument for which norm is returned.

Returns:

the L2 norm.