Comparison among criteria by standard test problems

• Solving two-dimensional problems in Standard Test Functions by the W4 method with the SV decomposition

Package Inclusion

• Sympy :: To define the set of nonlinear equations
• Numpy :: For iterative solver
• Scipy :: To use Linear solver
• Matplotlib :: For figure

import sympy as sy;
import numpy as np;
from sympy import exp,cos;
#from numpy import cos,sin,pi;
from scipy import linalg;
import matplotlib.pyplot as plt;

Parameter Setup

For the W4 method, we define

• dim : dimension of system (It determines the resolution)
• dtau : $\Delta\tau$ for the W4 evolution
• errmax : Criterion of error to stop iteration
• itermax : Maximum iteration.

dim=2;

dtau=0.8;
errmax = pow(10,-8);
itermax = pow(10,6);

### Variables
### x[1], x[2], x[3], ... , x[dim]
### 
x = []
for i in range(0,dim):
    vname = 'x[' + str(i) +']';
    x.append(sy.symbols(vname));

Definition of Source in the Sympy style

• F(x) : System of nonlinear equations
• |F(x)| : Absolute components
• J : Jacobian matrix $\partial$ F / $\partial$ x

f_rosenbrock = [];
f_rosenbrock.append( 10*(x[1]-x[0]**2) );
f_rosenbrock.append( 1 - x[0] );

f_freudenstein = [];
f_freudenstein.append( -13 +x[0] +((5-x[1])*x[1]-2)*x[1] );
f_freudenstein.append( -29 +x[0] +((x[1]+1)*x[1]-14)*x[1] );

f_powell = [];
f_powell.append( pow(10,4)*x[0]*x[1] -1 );
f_powell.append( exp(-x[0]) +exp(-x[1]) -1.0001 );

f_brown = [];
f_brown.append( x[0]*x[1]**2 -2*x[1] +x[0] -pow(10,6) );
f_brown.append( x[0]**2*x[1] -2*x[0] +x[1] -2*pow(10,-6) );

f_beale = [];
f_beale.append( 1.5 -x[0]*(1-x[1]) );
f_beale.append( 2.25 -x[0]*(1-x[1]**2) );

f_hueso = [];
f_hueso.append( (x[0]-1)**2*(x[0]-x[1]) );
f_hueso.append( (x[1]-2)**5*cos(2*x[0]/x[1]) );

f_fujisawa = [];
f_fujisawa.append( x[0]**2 +x[1]**2 -4 );
f_fujisawa.append( x[0]**2*x[1] -1 );

# input source
def calc_source(x,f):

### Summation of absolute values of all terms in each equation
###  to compare the error with the typical value of each equation
    fa = [];
    for j in range(dim):
        tmp = 0;
        for i in range(len(f[j].args)):
            tmp = tmp + abs(f[j].args[i]);
        fa.append(tmp);        

### Definition of variables, sources, absolute sources(x, F, |F|) as vector
    v = sy.Matrix([x]).transpose();
    F = sy.Matrix([f]).transpose();
    Fa = sy.Matrix([fa]).transpose();
    
### Analytic calculation of Jacobian(dim x dim Matrix)
    J = sy.Matrix();
    for i in range(dim):
        J = J.col_insert(i,F.diff(x[i]));


    fa1 = [];
    for j in range(dim):
        tmp = 0;
        for i in range(dim):
            tmp = tmp + abs(J[j,i]*v[i]);
        fa1.append(tmp);
    fa2 = [];
    for j in range(dim):
        tmp = abs(v[i]);
        fa2.append(tmp);

    Fa1 = sy.Matrix([fa1]).transpose();
    Fa2 = sy.Matrix([fa2]).transpose();
### From Sympy to Numpy
    arg = v.transpose();
    vn = sy.lambdify(arg, v, "numpy")
    Fn = sy.lambdify(arg, F, "numpy")
    Fan = sy.lambdify(arg, Fa, "numpy")
    Fa1n = sy.lambdify(arg, Fa1, "numpy")
    Fa2n = sy.lambdify(arg, Fa2, "numpy")
    Jacn = sy.lambdify(arg, J, "numpy")

    return vn,Fn,Fan,Fa1n,Fa2n,Jacn;

Transformation from Sympy style to Numpy style

#vn, Fn, Fan, Jacn = calc_source(x,f);
vn_rosenbrock, Fn_rosenbrock, Fan_rosenbrock, Fa1n_rosenbrock, Fa2n_rosenbrock, Jacn_rosenbrock = calc_source(x,f_rosenbrock);
vn_freudenstein, Fn_freudenstein, Fan_freudenstein, Fa1n_freudenstein, Fa2n_freudenstein, Jacn_freudenstein = calc_source(x,f_freudenstein);
vn_powell, Fn_powell, Fan_powell, Fa1n_powell, Fa2n_powell, Jacn_powell = calc_source(x,f_powell);
vn_brown, Fn_brown, Fan_brown, Fa1n_brown, Fa2n_brown, Jacn_brown = calc_source(x,f_brown);
vn_beale, Fn_beale, Fan_beale, Fa1n_beale, Fa2n_beale, Jacn_beale = calc_source(x,f_beale);
vn_hueso, Fn_hueso, Fan_hueso, Fa1n_hueso, Fa2n_hueso, Jacn_hueso = calc_source(x,f_hueso);
vn_fujisawa, Fn_fujisawa, Fan_fujisawa, Fa1n_fujisawa, Fa2n_fujisawa, Jacn_fujisawa = calc_source(x,f_fujisawa);

Definition :: W4SV method

def w4sv_analysis(v,F,Fa,Fa1,Fa2,J,ini,dt,itermax,errmax):
    
    p0 = np.zeros(dim);
    v0 = v(*ini).transpose()[0];

    evot = [];
#    evox = np.array([v0]).reshape(1,dim);
    F0=F(*v0).transpose()[0].transpose();
    evof = np.array([F0]).reshape(1,dim);
    F0a1=Fa(*v0).transpose()[0].transpose();
    F0a2=Fa1(*v0).transpose()[0].transpose();
    F0a3=Fa2(*v0).transpose()[0].transpose();
    evoa1 = np.array([F0a1]).reshape(1,dim);
    evoa2 = np.array([F0a2]).reshape(1,dim);
    evoa3 = np.array([F0a3]).reshape(1,dim);

    for i in range(itermax):
        J0=J(*v0);
        F0=F(*v0).transpose()[0].transpose();
        F0a1=Fa(*v0).transpose()[0].transpose();
        F0a2=Fa1(*v0).transpose()[0].transpose();
        F0a3=Fa2(*v0).transpose()[0].transpose();
#        print(v0);#,print(p0);
        err = 0.0;
        for k in range(dim):
            err = max(err,abs(F0[k]/F0a1[k]));

        if err < errmax:
            break;

        U, s, V = np.linalg.svd(J0);
        si = np.array([]);
        for j in range(dim):
            if s[j] > pow(10,-6):
                si = np.append(si,1.0/s[j]);
            else:
                si = np.append(si,1.0);

        S = np.diag(si);

        evot = np.append(evot,i);
#        evox = np.append(evox,np.array([v0]),axis=0);
#        evop = np.append(evop,np.array([p0]),axis=0);

        evof = np.append(evof,np.array([F0]),axis=0);
        evoa1 = np.append(evoa1,np.array([F0a1]),axis=0);
        evoa2 = np.append(evoa2,np.array([F0a2]),axis=0);
        evoa3 = np.append(evoa3,np.array([F0a3]),axis=0);

        srcx = V.transpose().dot(p0);
        tmp  = U.transpose().dot(F0);
        srcp = -2*p0 -S.dot(tmp);

        v0 = v0 +srcx*dt;
        p0 = p0 +srcp*dt;

    return v0, i, err, evot, evoa1, evoa2, evoa3, evof

Initial conditions for test problems

vini_rosenbrock = np.array([]);
vini_rosenbrock = np.append(vini_rosenbrock, 1.2);
vini_rosenbrock = np.append(vini_rosenbrock, 1.0);

vini_freudenstein = np.array([]);
vini_freudenstein = np.append(vini_freudenstein, 6.0);
vini_freudenstein = np.append(vini_freudenstein, 3.0);

vini_powell1 = np.array([]);
vini_powell1 = np.append(vini_powell1, 0.0);
vini_powell1 = np.append(vini_powell1, 1.0);

vini_powell2 = np.array([]);
vini_powell2 = np.append(vini_powell2, 1.0);
vini_powell2 = np.append(vini_powell2, 1.0);


vini_brown = np.array([]);
vini_brown = np.append(vini_brown, 1.0);
vini_brown = np.append(vini_brown, 1.0);

vini_beale1 = np.array([]);
vini_beale1 = np.append(vini_beale1, 1.0);
vini_beale1 = np.append(vini_beale1, 1.0);

vini_beale2 = np.array([]);
vini_beale2 = np.append(vini_beale2, 0.0);
vini_beale2 = np.append(vini_beale2, 2.0);


vini_hueso = np.array([]);
vini_hueso = np.append(vini_hueso, 1.5);
vini_hueso = np.append(vini_hueso, 2.5);

vini_fujisawa1 = np.array([]);
vini_fujisawa1 = np.append(vini_fujisawa1, 0.0);
vini_fujisawa1 = np.append(vini_fujisawa1, 1.0);

vini_fujisawa2 = np.array([]);
vini_fujisawa2 = np.append(vini_fujisawa2, 0.0);
vini_fujisawa2 = np.append(vini_fujisawa2, -1.0);

Demonstrations

• Rosenbrock's problem
• Freudenstein's problem
• Powell's problem (Two initial conditions)
• Brown's problem
• Beale's problem (Two initial conditions)
• Hueso's problem
• Fujisawa's problem (Two initial conditions)

Outputs :

( Solution x and y, Number of iterations to obtain the solution, Convergence Error )

vans_rosenbrock, iter_rosenbrock, err, t0_rosenbrock, a1_rosenbrock, a2_rosenbrock, a3_rosenbrock, a0_rosenbrock = w4sv_analysis(vn_rosenbrock,Fn_rosenbrock,Fan_rosenbrock,Fa1n_rosenbrock,Fa2n_rosenbrock,Jacn_rosenbrock,vini_rosenbrock,dtau,itermax,errmax);
vans_rosenbrock, iter_rosenbrock, err

(array([1.        , 0.99999999]), 14, 8.874838939484052e-09)

vans_freudenstein, iter_freudenstein, err, t0_freudenstein, a1_freudenstein, a2_freudenstein, a3_freudenstein, a0_freudenstein = w4sv_analysis(vn_freudenstein,Fn_freudenstein,Fan_freudenstein, Fa1n_freudenstein, Fa2n_freudenstein,Jacn_freudenstein,vini_freudenstein,dtau,itermax,errmax);
vans_freudenstein, iter_freudenstein, err

(array([4.99999995, 4.00000001]), 18, 3.671612156005402e-09)

vans_powell1, iter_powell1, err, t0_powell1, a1_powell1, a2_powell1, a3_powell1, a0_powell1 = w4sv_analysis(vn_powell,Fn_powell,Fan_powell,Fa1n_powell,Fa2n_powell,Jacn_powell,vini_powell1,dtau,itermax,errmax);
vans_powell1, iter_powell1, err

(array([1.09815945e-05, 9.10614562e+00]), 33, 5.128652417964679e-09)

vans_powell2, iter_powell2, err, t0_powell2, a1_powell2, a2_powell2, a3_powell2, a0_powell2 = w4sv_analysis(vn_powell,Fn_powell,Fan_powell,Fa1n_powell,Fa2n_powell,Jacn_powell,vini_powell2,dtau,itermax,errmax);
vans_powell2, iter_powell2, err

(array([1.09815952e-05, 9.10614497e+00]), 66, 8.089358248742966e-09)

vans_brown, iter_brown, err, t0_brown, a1_brown, a2_brown, a3_brown, a0_brown = w4sv_analysis(vn_brown,Fn_brown,Fan_brown,Fa1n_brown,Fa2n_brown,Jacn_brown,vini_brown,dtau,itermax,errmax);
vans_brown, iter_brown, err

(array([1.00000000e+06, 2.00000002e-06]), 10724, 6.0799684975196886e-09)

vans_beale1, iter_beale1, err, t0_beale1, a1_beale1, a2_beale1, a3_beale1, a0_beale1 = w4sv_analysis(vn_beale,Fn_beale,Fan_beale,Fa1n_beale,Fa2n_beale,Jacn_beale,vini_beale1,dtau,itermax,errmax);
vans_beale1, iter_beale1, err

(array([2.99999995, 0.49999999]), 29, 3.878248453483837e-09)

vans_beale2, iter_beale2, err, t0_beale2, a1_beale2, a2_beale2, a3_beale2, a0_beale2 = w4sv_analysis(vn_beale,Fn_beale,Fan_beale,Fa1n_beale,Fa2n_beale,Jacn_beale,vini_beale2,dtau,itermax,errmax);
vans_beale2, iter_beale2, err

(array([2.99999995, 0.49999999]), 32, 3.934620447417195e-09)

vans_hueso, iter_hueso, err, t0_hueso, a1_hueso, a2_hueso, a3_hueso, a0_hueso = w4sv_analysis(vn_hueso,Fn_hueso,Fan_hueso,Fa1n_hueso,Fa2n_hueso,Jacn_hueso,vini_hueso,dtau,itermax,errmax);
vans_hueso, iter_hueso, err

(array([1.00006356, 2.02433785]), 34, 8.539072485799175e-09)

vans_fujisawa1, iter_fujisawa1, err, t0_fujisawa1, a1_fujisawa1, a2_fujisawa1, a3_fujisawa1, a0_fujisawa1 = w4sv_analysis(vn_fujisawa,Fn_fujisawa,Fan_fujisawa,Fa1n_fujisawa,Fa2n_fujisawa,Jacn_fujisawa,vini_fujisawa1,dtau,itermax,errmax);
vans_fujisawa1, iter_fujisawa1, err

(array([-0.73307679,  1.86080586]), 16, 3.8859071365125804e-09)

vans_fujisawa2, iter_fujisawa2, err, t0_fujisawa2, a1_fujisawa2, a2_fujisawa2, a3_fujisawa2, a0_fujisawa2 = w4sv_analysis(vn_fujisawa,Fn_fujisawa,Fan_fujisawa,Fa1n_fujisawa,Fa2n_fujisawa,Jacn_fujisawa,vini_fujisawa2,dtau,itermax,errmax);
vans_fujisawa2, iter_fujisawa2, err

(array([-1.98379242,  0.25410168]), 82, 6.008461470011285e-09)

Comparison of criteria by Fujisawa's Problem

• Newton-Raphson method cannot find the solution starting from this initial guess

f1a0 = [];
for i in range(0,iter_fujisawa2):
    f1a0 = np.append(f1a0,a0_fujisawa2[i][0]);

f1a1 = [];
for i in range(0,iter_fujisawa2):
    f1a1 = np.append(f1a1,a1_fujisawa2[i][0]);
    
f1a2 = [];
for i in range(0,iter_fujisawa2):
    f1a2 = np.append(f1a2,a2_fujisawa2[i][0]);

f1a3 = [];
for i in range(0,iter_fujisawa2):
    f1a3 = np.append(f1a3,a3_fujisawa2[i][0]);

f2a0 = [];
for i in range(0,iter_fujisawa2):
    f2a0 = np.append(f2a0,a0_fujisawa2[i][1]);

f2a1 = [];
for i in range(0,iter_fujisawa2):
    f2a1 = np.append(f2a1,a1_fujisawa2[i][1]);
    
f2a2 = [];
for i in range(0,iter_fujisawa2):
    f2a2 = np.append(f2a2,a2_fujisawa2[i][1]);

f2a3 = [];
for i in range(0,iter_fujisawa2):
    f2a3 = np.append(f2a3,a3_fujisawa2[i][1]);
    
t = t0_fujisawa2

### Figures
fig = plt.figure(figsize=(12,3));
ax1 = fig.add_subplot(1,2,1);
ax1.plot(t,(abs(f1a0)/f1a1),marker='o',label='Definition 1');
ax1.plot(t,(abs(f1a0)/f1a2),marker='x',label='Definition 2');
ax1.plot(t,(abs(f1a0)/f1a3),marker='d',label='Definition 3');
ax1.set_xlabel("$Iteration$");
ax1.set_ylabel("$|F_1|/|F_1|_{crit}$");
ax1.set_yscale('log');
ax1.set_title('Error evolution of $F_1$')
ax1.legend();

ax2 = fig.add_subplot(1,2,2);
ax2.plot(t,(abs(f2a0)/f2a1),marker='o',label='Definition 1');
ax2.plot(t,(abs(f2a0)/f2a2),marker='x',label='Definition 2');
ax2.plot(t,(abs(f2a0)/f2a3),marker='d',label='Definition 3');
ax2.set_xlabel("$Iteration$");
ax2.set_ylabel("$|F_2|/|F_2|_{crit}$");
ax2.set_yscale('log');
ax2.set_title('Error evolution of $F_{2}$')
ax2.legend();

fig.tight_layout();

<ipython-input-22-764b841bd1b9>:15: RuntimeWarning: divide by zero encountered in true_divide
  ax2.plot(t,(abs(f2a0)/f2a2),marker='x',label='Definition 2');

Comparison of criteria by Hueso's Problem

• The Jacobian matrix for this problem at the solution is singular.

f1a0 = [];
for i in range(0,iter_hueso):
    f1a0 = np.append(f1a0,a0_hueso[i][0]);

f1a1 = [];
for i in range(0,iter_hueso):
    f1a1 = np.append(f1a1,a1_hueso[i][0]);
    
f1a2 = [];
for i in range(0,iter_hueso):
    f1a2 = np.append(f1a2,a2_hueso[i][0]);

f1a3 = [];
for i in range(0,iter_hueso):
    f1a3 = np.append(f1a3,a3_hueso[i][0]);

f2a0 = [];
for i in range(0,iter_hueso):
    f2a0 = np.append(f2a0,a0_hueso[i][1]);

f2a1 = [];
for i in range(0,iter_hueso):
    f2a1 = np.append(f2a1,a1_hueso[i][1]);
    
f2a2 = [];
for i in range(0,iter_hueso):
    f2a2 = np.append(f2a2,a2_hueso[i][1]);

f2a3 = [];
for i in range(0,iter_hueso):
    f2a3 = np.append(f2a3,a3_hueso[i][1]);
    
t = t0_hueso

### Figures
fig = plt.figure(figsize=(12,3));
ax1 = fig.add_subplot(1,2,1);
ax1.plot(t,(abs(f1a0)/f1a1),marker='o',label='Definition 1');
ax1.plot(t,(abs(f1a0)/f1a2),marker='x',label='Definition 2');
ax1.plot(t,(abs(f1a0)/f1a3),marker='d',label='Definition 3');
ax1.set_xlabel("$Iteration$");
ax1.set_ylabel("$|F_1|/|F_1|_{crit}$");
ax1.set_yscale('log');
ax1.set_title('Error evolution of $F_1$')
ax1.legend();

ax2 = fig.add_subplot(1,2,2);
ax2.plot(t,(abs(f2a0)/f2a1),marker='o',label='Definition 1');
ax2.plot(t,(abs(f2a0)/f2a2),marker='x',label='Definition 2');
ax2.plot(t,(abs(f2a0)/f2a3),marker='d',label='Definition 3');
ax2.set_xlabel("$Iteration$");
ax2.set_ylabel("$|F_2|/|F_2|_{crit}$");
ax2.set_yscale('log');
ax2.set_title('Error evolution of $F_{2}$')
ax2.legend();

fig.tight_layout();

Comments on the first Figures: Fujisawa's problem

• $|F_i|_{crit}$ is defined in the three ways:
  1. the sum of absolute values of all terms,
  2. the sum of Jacobian times variable, $J_{ij}x_{j}$
  3. the absolute value of $x_{i}$
• For Fujisawa's problem, any definition of the criterion to stop the W4 iteration works well.
• We sometimes face with the singularity that makes the second criterion worse to evaluate the order of equation.

Comments on the second Figures: Hueso's problem

• For Hueso's problem, the second definition becomes worse as we reach the solution.
• To conclude, it seems practically fine to evaluate the order of equations by the third definition, which is just to use the order of variables.
• Note that the third definition may work unless the solution is exactly zero.