Comparison among different criteria to stop W4 iteration

• A criterion to stop the iteration is necessary.
• From the point of view of implementation, possible criteria will be compared.

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 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 evolution
• errmax : Criterion of error to stop iteration
• itermax : Maximum iteration.

dim=10;

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

### 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

$\frac{d^2\theta}{d\xi^2}+\frac{2}{\xi}\frac{d\theta}{d\xi}+r_{*}^2\theta^{n}=0$

### Source
### f[0], f[1], f[2], ... , f[dim-1]
n = 1;
f = [];

### Coordinate xi
dx = 1.0/dim;
xi = [];
for i in range(dim+1):
    xi.append( dx*i );

    
### Source for Eigenvalue x[0]
f.append( 6*(x[1] -1)/dx**2  +x[0] );

### Source for x[1] (theta(dx)) with boundary condition theta(0)=1
f.append( (x[2] -2*x[1] +1)/dx**2 +(x[2]-1)/xi[1]/dx +x[1]**n*x[0] );

for i in range(2,dim-1):
    f.append( (x[i+1] -2*x[i] +x[i-1])/dx**2 +(x[i+1]-x[i-1])/xi[i]/dx +x[i]**n*x[0] );

### Source for x[dim-1] (theta(1-dx)) with boundary condition theta(1)=0
f.append( (0 -2*x[dim-1] +x[dim-2])/dx**2 +(0-x[dim-2])/xi[dim-1]/dx +x[dim-1]**n*x[0] );

#f
tmp = f[0];
#norm = 0.0;
for i in range(1,dim):
    tmp = tmp +f[i]*x[i];
#    norm = norm +x[i]**n*x[i];

f[0] = tmp;

# 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
    fa1 = [];
    for j in range(dim):
        tmp = 0;
        for i in range(len(f[j].args)):
            tmp = tmp + abs(f[j].args[i]);
        fa1.append(tmp);        

### Definition of variables, sources, absolute sources(x, F, |F|) as vector
    v = sy.Matrix([x]).transpose();
    F = sy.Matrix([f]).transpose();
    
    fa3 = [];
    for j in range(dim):
        tmp = abs(v[j]);
        fa3.append(tmp);        
        
### 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]));    

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

    Fa1 = sy.Matrix([fa1]).transpose();
    Fa2 = sy.Matrix([fa2]).transpose();
    Fa3 = sy.Matrix([fa3]).transpose();

### From Sympy to Numpy
    arg = v.transpose();
    vn = sy.lambdify(arg, v, "numpy")
    Fn = sy.lambdify(arg, F, "numpy")
    Fa1n = sy.lambdify(arg, Fa1, "numpy")
    Fa2n = sy.lambdify(arg, Fa2, "numpy")
    Fa3n = sy.lambdify(arg, Fa3, "numpy")
    Jacn = sy.lambdify(arg, J, "numpy")        
        
    return vn,Fn,Fa1n,Jacn,Fa2n,Fa3n;

vn, Fn, Fa1n, Jacn, Fa2n, Fa3n = calc_source(x,f);

Definition :: W4LH method

### Analysis by the W4 method with the LH decomposition
def w4analysis2(v,F,Fa1,Fa2,Fa3,J,ini,dt,itermax,errmax):
    
### Initialization for x and p as v0 and p0
    p0 = np.zeros(dim);
    v0 = v(*ini).transpose()[0];

    evot = [];
#    evot = np.append(evot,0);
    evox = np.array([v0]).reshape(1,dim);
    evop = np.array([p0]).reshape(1,dim);
    evolam = [];
    J0=J(*v0);
    lam = linalg.eigvals(J0)
    lmin=abs(lam[0]);
    lmax=abs(lam[0]);
    for d in range(dim):
        if lmin > abs(lam[d]):
            lmin = abs(lam[d]);
        if lmax < abs(lam[d]):
            lmax = abs(lam[d]);

    evolam = np.append(evolam,lmax/lmin);

    
    F0=F(*v0).transpose()[0].transpose();
    evof = np.array([F0]).reshape(1,dim);
    F0a1=Fa1(*v0).transpose()[0].transpose();
    F0a2=Fa2(*v0).transpose()[0].transpose();
    F0a3=Fa3(*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);
    
    
### Main Iteration Loop
    for i in range(itermax):
        J0=J(*v0);
        F0=F(*v0).transpose()[0].transpose();
        F0a1=Fa1(*v0).transpose()[0].transpose();
        F0a2=Fa2(*v0).transpose()[0].transpose();
        F0a3=Fa3(*v0).transpose()[0].transpose();
        
### Criterion to stop the iteration
### Evaluation of Source F
        err = 0.0;
        for k in range(dim):
            err = max(err,abs(F0[k]/F0a1[k]));
            
        if err < errmax:
#            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);
            break;

### Decomposition of Jacobian into L and H
        Q, R = linalg.qr(J0.transpose());
        srcp = -2*p0 -linalg.solve(R.transpose(),F0);
        srcx = Q.dot(p0);
        v0 = v0 +srcx*dt;
        p0 = p0 +srcp*dt;
        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);

        lam = linalg.eigvals(J0)
        lmin=abs(lam[0]);
        lmax=abs(lam[0]);
        for d in range(dim):
            if lmin > abs(lam[d]):
                lmin = abs(lam[d]);
            if lmax < abs(lam[d]):
                lmax = abs(lam[d]);

        evolam = np.append(evolam,lmax/lmin);

    return v0, i, err, evolam, evot, evox, evop, evoa1, evoa2, evoa3, evof

### Initial guess
vini = np.array([]);
vini = np.append(vini, 1.);
for i in range(1,dim):
    vini = np.append(vini, (1.0 -dx*i));

### Nonlinear Solver by the W4 method
vans, iter, err, lam, t, xr, pr, a1, a2, a3, a0 = w4analysis2(vn,Fn,Fa1n,Fa2n,Fa3n,Jacn,vini,dtau,itermax,errmax);

### Basic information
print("Dimension : ",dim);
print("Initial Condition number : ",lam[1]);
print("Iteration : ",iter,", error : ", err);
if iter < itermax-1:
    print("< Solved. >");
else:
    print("< Failed. >");

Dimension :  10
Initial Condition number :  20.20479855607713
Iteration :  15 , error :  5.511184118127499e-09
< Solved. >

### Solution
sol = [];
sol = np.append(sol,1.0);
for i in range(1,dim):
    sol = np.append(sol,xr[iter][i]);
    
sol = np.append(sol,0.0);

r = [];
for i in range(dim+1):
    r = np.append(r,xi[i]*xr[iter][0]**0.5);
    
### Initial guess    
sola = [];
sola = np.append(sola,1.0);
for i in range(1,dim):
    sola = np.append(sola,vini[i]);
    
sola = np.append(sola,0.0);

ra = [];
for i in range(dim+1):
    ra = np.append(ra,xi[i]*vini[0]**0.5);

### Figures
fig = plt.figure(figsize=(12,5));
ax1 = fig.add_subplot(1,1,1);
ax1.plot(ra,sola,marker='x',label='Initial guess');
ax1.plot(r,sol,marker='o',label='Solution');
ax1.set_xlabel("$r$");
ax1.set_ylabel("Solution $\Theta$");

ax1.legend();
fig.tight_layout();

Comments on the setup to solve Lane-Emden Equation

What is actually solved?

• To solve the differential equation, we discretize it with the number of dimension, which the space is divided.
• Variables are solved using discretized equations starting from the initial guess.

f1a0 = [];
for i in range(0,iter):
    f1a0 = np.append(f1a0,a0[i][1]);

f1a1 = [];
for i in range(0,iter):
    f1a1 = np.append(f1a1,a1[i][1]);
    
f1a2 = [];
for i in range(0,iter):
    f1a2 = np.append(f1a2,a2[i][1]);

f1a3 = [];
for i in range(0,iter):
    f1a3 = np.append(f1a3,a3[i][1]);

f2a0 = [];
for i in range(0,iter):
    f2a0 = np.append(f2a0,a0[i][dim//2]);
    
f2a1 = [];
for i in range(0,iter):
    f2a1 = np.append(f2a1,a1[i][dim//2]);
    
f2a2 = [];
for i in range(0,iter):
    f2a2 = np.append(f2a2,a2[i][dim//2]);

f2a3 = [];
for i in range(0,iter):
    f2a3 = np.append(f2a3,a3[i][dim//2]);

f3a0 = [];
for i in range(0,iter):
    f3a0 = np.append(f3a0,a0[i][dim-1]);
    
f3a1 = [];
for i in range(0,iter):
    f3a1 = np.append(f3a1,a1[i][dim-1]);
    
f3a2 = [];
for i in range(0,iter):
    f3a2 = np.append(f3a2,a2[i][dim-1]);

f3a3 = [];
for i in range(0,iter):
    f3a3 = np.append(f3a3,a3[i][dim-1]);

### Figures
fig = plt.figure(figsize=(12,3));
ax1 = fig.add_subplot(1,3,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 at $\Theta_{1}$')
ax1.legend();

ax2 = fig.add_subplot(1,3,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_{N/2}|/|F_{N/2}|_{crit}$");
ax2.set_yscale('log');
ax2.set_title('Error evolution at $\Theta_{N/2}$')
ax2.legend();

ax3 = fig.add_subplot(1,3,3);
ax3.plot(t,(abs(f3a0)/f3a1),marker='o',label='Definition 1');
ax3.plot(t,(abs(f3a0)/f3a2),marker='x',label='Definition 2');
ax3.plot(t,(abs(f3a0)/f3a3),marker='d',label='Definition 3');
ax3.set_xlabel("$Iteration$");
ax3.set_ylabel("$|F_{N}|/|F_{N}|_{crit}$");
ax3.set_yscale('log');
ax3.set_title('Error evolution at $\Theta_{N}$')
ax3.legend();

fig.tight_layout();

Comments on the criteria of W4 iteration for LE equation

How to see the figures

• Each figure shows the evolution of error defined as $|F_i|/|F_i|_{crit}$, which is the absolute value of the discretized equation divided by a typical value of equation.
• $i$ denotes the number of grid points.
• $|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}$
• From this example, the definition $2$ seems working well, but it is noted that the definition $3$ is not so bad.