[SciPy-User] solving integro-differential equation over a square with variable boundary condition on one edge (scipy.optimize)
Edrisse Chermak
edrisse.chermak at gmail.com
Wed Jan 21 00:29:09 EST 2015
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>
More information about the SciPy-User
mailing list