Nonlinear Solver Lab
Demonstration of the W4SV method
- To demonstrate the W4 method, we show below how to solve the system of nonlinear equations by the Python version of the W4 method with the SV decomposition.
- In this example, standard test problems in the paper are solved. You can obtain the Python notebook in https://hir0ok.github.io/w4/W4demo_W4SV_TestFunctions.ipynb .
W4SV demonstration¶
- 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
In [2]:
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.
In [3]:
dim=2;
In [4]:
dtau=0.8;
errmax = pow(10,-4);
itermax = pow(10,7);
In [5]:
### 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
In [6]:
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 );
In [7]:
# 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]));
### 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")
Jacn = sy.lambdify(arg, J, "numpy")
return vn,Fn,Fan,Jacn;
Transformation from Sympy style to Numpy style¶
In [8]:
#vn, Fn, Fan, Jacn = calc_source(x,f);
vn_rosenbrock, Fn_rosenbrock, Fan_rosenbrock, Jacn_rosenbrock = calc_source(x,f_rosenbrock);
vn_freudenstein, Fn_freudenstein, Fan_freudenstein, Jacn_freudenstein = calc_source(x,f_freudenstein);
vn_powell, Fn_powell, Fan_powell, Jacn_powell = calc_source(x,f_powell);
vn_brown, Fn_brown, Fan_brown, Jacn_brown = calc_source(x,f_brown);
vn_beale, Fn_beale, Fan_beale, Jacn_beale = calc_source(x,f_beale);
vn_hueso, Fn_hueso, Fan_hueso, Jacn_hueso = calc_source(x,f_hueso);
vn_fujisawa, Fn_fujisawa, Fan_fujisawa, Jacn_fujisawa = calc_source(x,f_fujisawa);
Definition :: W4SV method¶
In [9]:
def w4sv(v,F,Fa,J,ini,dt,itermax,errmax):
p0 = np.zeros(dim);
v0 = v(*ini).transpose()[0];
for i in range(itermax):
J0=J(*v0);
F0=F(*v0).transpose()[0].transpose();
Fa0=Fa(*v0).transpose()[0].transpose();
# print(v0);#,print(p0);
err = 0.0;
for k in range(dim):
err = max(err,abs(F0[k]/Fa0[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);
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
Initial conditions for test problems¶
In [10]:
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 )¶
In [11]:
vans_rosenbrock, iter, err = w4sv(vn_rosenbrock,Fn_rosenbrock,Fan_rosenbrock,Jacn_rosenbrock,vini_rosenbrock,dtau,itermax,errmax);
vans_rosenbrock, iter, err
Out[11]:
In [12]:
vans_freudenstein, iter, err = w4sv(vn_freudenstein,Fn_freudenstein,Fan_freudenstein,Jacn_freudenstein,vini_freudenstein,dtau,itermax,errmax);
vans_freudenstein, iter, err
Out[12]:
In [13]:
vans_powell1, iter, err = w4sv(vn_powell,Fn_powell,Fan_powell,Jacn_powell,vini_powell1,dtau,itermax,errmax);
vans_powell1, iter, err
Out[13]:
In [21]:
vans_powell2, iter, err = w4sv(vn_powell,Fn_powell,Fan_powell,Jacn_powell,vini_powell2,dtau,itermax,errmax);
vans_powell2, iter, err
Out[21]:
In [22]:
vans_brown, iter, err = w4sv(vn_brown,Fn_brown,Fan_brown,Jacn_brown,vini_brown,dtau,itermax,errmax);
vans_brown, iter, err
Out[22]:
In [23]:
vans_beale1, iter, err = w4sv(vn_beale,Fn_beale,Fan_beale,Jacn_beale,vini_beale1,dtau,itermax,errmax);
vans_beale1, iter, err
Out[23]:
In [24]:
vans_beale2, iter, err = w4sv(vn_beale,Fn_beale,Fan_beale,Jacn_beale,vini_beale2,dtau,itermax,errmax);
vans_beale2, iter, err
Out[24]:
In [25]:
vans_hueso, iter, err = w4sv(vn_hueso,Fn_hueso,Fan_hueso,Jacn_hueso,vini_hueso,dtau,itermax,errmax);
vans_hueso, iter, err
Out[25]:
In [19]:
vans_fujisawa1, iter, err = w4sv(vn_fujisawa,Fn_fujisawa,Fan_fujisawa,Jacn_fujisawa,vini_fujisawa1,dtau,itermax,errmax);
vans_fujisawa1, iter, err
Out[19]:
In [26]:
vans_fujisawa2, iter, err = w4sv(vn_fujisawa,Fn_fujisawa,Fan_fujisawa,Jacn_fujisawa,vini_fujisawa2,dtau,itermax,errmax);
vans_fujisawa2, iter, err
Out[26]:
In [ ]: