[SciPy-User] solving integro-differential equation over a square with variable boundary condition on one edge (scipy.optimize)

[Edrisse Chermak]

Dear Scipy users,

I'm trying to solve a nonlinear integrodifferential equation over a 
square with Scipy by introducing a *variable* boundary condition on one 
edge of the square.Scipy documentation example on nonlinear solvers 
<http://docs.scipy.org/doc/scipy-0.14.0/reference/optimize.nonlin.html> 
[1] gives a sample code on how to do solve this equation, but with a set 
of 4 *constant* boundary conditions on each edge : P(x,1) = 1 and P=0 on 
any other boundary of the square [0,1]x[0,1]. That means the top edge 
boundary condition is P=1 and P=0 for the 3 other edges :

|from  scipy.optimizeimport  newton_krylov
from  numpyimport  cosh,  zeros_like,  mgrid,  zeros

# Defining 2D grid of 75x75 dots :
nx,  ny=  75,  75
# Defining grid steps
hx,  hy=  1./(nx-1),  1./(ny-1)
# Boundary Conditions for function P(x,y) over square [0,1]x[0,1]
P_left,  P_right=  0,  0                         # P(1,y)=0 and P(0,y)=0
P_top,  P_bottom=  1,  0                         # P(x,1)=1 and P(x,0)=0

# Defining integro-differential equation for function P(x,y)
def  residual(P):
     d2x=  zeros_like(P)
     d2y=  zeros_like(P)
     # Central Difference Approximation
     d2x[1:-1]  =  (P[2:]    -  2*P[1:-1]  +  P[:-2])  /  hx/hx
     d2x[0]     =  (P[1]     -  2*P[0]     +  P_left)/hx/hx
     d2x[-1]    =  (P_right-  2*P[-1]    +  P[-2])/hx/hx
     d2y[:,1:-1]  =  (P[:,2:]  -  2*P[:,1:-1]  +  P[:,:-2])/hy/hy
     d2y[:,0]     =  (P[:,1]   -  2*P[:,0]     +  P_bottom)/hy/hy
     d2y[:,-1]    =  (P_top-  2*P[:,-1]    +  P[:,-2])/hy/hy
     # expression of the integro-differential equation
     return  d2x+  d2y-  10*cosh(P).mean()**2

# Defining a guess starting solution
guess=  zeros((nx,  ny),  float)
# Iterating to find the solution
sol=  newton_krylov(residual,  guess,  method='lgmres',  verbose=1)
print('Residual: %g'  %  abs(residual(sol)).max())
# visualize the solution
import  matplotlib.pyplotas  plt
x,  y=  mgrid[0:1:(nx*1j),  0:1:(ny*1j)]
plt.pcolor(x,  y,  sol)
plt.colorbar()
plt.show()|

I want to solve the same problem, with a *variable* boundary condition 
on the top edge that is:

|P_top=  P(x,1)  =  1/(0.5  -  abs(x))|

That means that P would have a different value on the top edge of the 
square, as shown on the following figure :

         y
           /\ *(P_top)*
           | *P(x,1)=1/(0.5-|x|)*
          1|__________________
           |                  |
           |                  |
P(0,y)=0  |                  | P(1,y)=0
(P_left)  |                  | (P_right)
           |                  |
           |                  |
           |__________________|_____> x
          0                   1
                P(x,0)=0
               (P_bottom)

I tried to insert P_top with the formula shown on line 11 of below code 
but I got an error message that variable 'x' is not defined. Would 
someone have any idea on how to define properly this P_top variable 
boundary condition in the residual ?

Thanks in advance,

[1] http://docs.scipy.org/doc/scipy-0.14.0/reference/optimize.nonlin.html

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20150121/398fd80c/attachment.html>