0

I am solving for f(kx,ky) that satisfies following integral equation, k is vector in 2d

equation:

here is the equation

I am mapping f(kx,ky) into 2d N*N grid and solving a system of total N^2 equations.

from scipy.optimize import newton_krylov
import numpy as n
L=10.0     # cut off of integral
N=100       # number of devision from 0 to L
delta=L/(N-1)

def f0(f):                # define a larger matrix used in summation
    f1=n.ones((2*N,2*N))
    for i in range(0, N):
       for j in range(0, N):
           f1[i][j]=f[i][j]
    return f1       


def equation(f):
    f2=n.empty((N,N))
    for i in range(0, N):      # i,j index of k
       for j in range(0, N):
          sum=0
          for a in range(-N+1, N):       # a,b index of p
             for b in range(-N+1, N):
                if a!=0 or b!=0:
                   for c in range(-N+1, N):       # c,d index of q                          
                      for d in range(-N+1, N):
                         if c!=0 or d!=0:
                            if a!=c or b!=d:
                               if c!=-i or d!=-j:
                                  sum=sum-delta**2/(delta**2*(a**2+b**2)\
 +f0(f)[abs(a)][abs(b)])/(delta**2*((a-c)**2+(b-d)**2)\
 +f0(f)[abs(a-c)][abs(b-d)])*(1/(delta**2*((c+i)**2+(d+j)**2)\
 +f0(f)[abs(c+j)][abs(d+j)])-1/(delta**2*(c**2+d**2)\
 +f0(f)[abs(c)][abs(d)]))
          f2[i][j]=sum-f[i][j]
    return f2       

guess =n.ones((N, N), float)
sol = newton_krylov(equation, guess, method='lgmres', verbose=1)

print(sol)

I first tried with very small parameter l=1, N=4. The solution finding is very slow. 

0:  |F(x)| = 23.6576; step 1; tol 0.487986
.....
519:  |F(x)| = 21.4812; step 3.15471e-09; tol 0.9
.....
957:  |F(x)| = 20.4225; step 2.03396e-06; tol 0.899997
.....

and still running.

The initial guess is far from the actual solution maybe that is why. Still I want to know whether it is a good way to use newton_krylov in my case?

Also noticing f(0) should be equal to zero in the actual solution. Is it reasonable to do the summation with exclusion of singular point as the integral?

One more question, how to force f(0,0) in the solution always to be equal to zero?

B. Desai
  • 16,414
  • 5
  • 26
  • 47
p.s
  • 51
  • 3
  • 2
    Four instances of `f0(f)[X][Y]` in the innermost loop are _really_ killing you. You create four of them at each iteration and then immediately dispose of them. Calculate `g=f0(f)` before entering the outermost loop and use it instead of `f0(f)` throughout the function `equation`. – DYZ Sep 21 '17 at 04:52
  • I could have use f directly for equation(f). However f is defined as N*N matrix and equation(f) needs the value of it a larger domain. So I defined a new matrix has same value as f except it is larger. There must be smarter ways to do it. – p.s Sep 21 '17 at 05:41
  • It's ok to define a bigger matrix. Just don't do it four time per innermost iteration. – DYZ Sep 21 '17 at 05:46
  • this sort of stuff is really in the domain of Cython and Numba – Eli Korvigo Sep 21 '17 at 06:03

0 Answers0