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

[Hanno Klemm]

, please excuse my brevity. 

> On 21.01.2015, at 06:29, Edrisse Chermak <edrisse.chermak at gmail.com> wrote:
> 
> 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 [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.optimize import newton_krylov
> from numpy import 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.pyplot as 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
> 
> _______________________________________________
> SciPy-User mailing list
> SciPy-User at scipy.org
> http://mail.scipy.org/mailman/listinfo/scipy-user

P_top =1./(0.5 - np.linspace(0,1,nx))

should do the trick. For the interval 0,1 you don't need to take the modulus. 

Best,
Hanno

hanno.klemm at me.com

Sent from my mobile device
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20150121/d07a2a06/attachment.html>