"""
Solving:
y'' = y'^2
y(0.5) = 0
y(1) = -4
Exact solution:
-ln(2x*exp(4)-exp(4)-2x+2)
Variational Form:
a(u',v') = -int( u'* v')
l(v) = int( u'**2 * v)
"""
from dolfin import *
import numpy
mesh = IntervalMesh(20, 0.5, 1)
V = FunctionSpace(mesh, 'CG', 2)
# defining boundary conditions
ul =Expression("0.5-x[0]")
def Boundary(x, on_boundary):
return on_boundary
bc1 = DirichletBC(V, ul, Boundary)
#defining the right boundary
ur = Constant(-4)
def Boundary(x, on_boundary):
return on_boundary
bc2 = DirichletBC(V, ur, Boundary)
bcs = [bc1, bc2]
# Setting up variatinal Problem
u = TrialFunction(V)
v =TestFunction(V)
u_k = interpolate(Constant(1.0), V)
fd = Dx(u_k, 0)
a = -inner(nabla_grad(u), nabla_grad(v))*dx
L = v*fd**2*dx
# setting up picard iteration
u = Function(V)
tol = 1.0E-2
eps = 1.0
maxiter = 10000
iteration = 0
while iteration < maxiter:
iteration += 1
solve(a ==L, u, bcs)
diff = u.vector().array() - u_k.vector().array()
eps = numpy.linalg.norm(diff, ord = numpy.Inf)
if eps < tol:
print "Procedure successful"
plot(u, interactive = True)
quit()
else:
u_k.assign(u)
print "Max. number of iterations have reached."
What am I doing wrong?
