import numpy as np
from matplotlib import pyplot as plt
def yexact(x):
return 1.0/(np.exp(x)+np.exp(-x))
def f(y,dy):
return -y + 2.0*dy**2/y
def phi(y):
n = len(y) - 1
h = 2.0/n
res = np.zeros(n-1)
for i in range(1,n):
dy = (y[i+1] - y[i-1])/(2.0*h)
res[i-1] = (y[i-1]-2.0*y[i]+y[i+1])/h**2 - f(y[i],dy)
return res
def dphi(y):
n = len(y) - 1
h = 2.0/n
res = np.zeros((n-1,n-1))
for i in range(1,n):
dy = (y[i+1] - y[i-1])/(2.0*h)
if i > 1:
res[i-1,i-2] = 1.0/h**2 + 4.0*dy/y[i]
res[i-1,i-1] =-2.0/h**2 + 1.0 - 2.0*(dy/y[i])**2
if i < n-1:
res[i-1,i ] = 1.0/h**2 - 4.0*dy/y[i]
return res
n = 100
# Initial guess for y
y = np.zeros(n+1)
y[:] = 1.0/(exp(1) + exp(-1))
maxiter = 100
TOL = 1.0e-6
it = 0
while it < maxiter:
b = phi(y)
if np.linalg.norm(b) < TOL:
break
A = dphi(y)
v = np.linalg.solve(A,b)
y[1:n] = y[1:n] - v
it += 1
print "Number of iterations = %d" % it
x = np.linspace(-1.0,1.0,n+1)
plt.plot(x,y,'o',x,yexact(x))
plt.legend(("Numerical","Exact"));
Number of iterations = 13