In [1]:
import sympy as sp
sp.init_printing()


## Make up the BVP¶

Choose the LHS and the solution and compute the RHS.

In [2]:
x = sp.symbols('x')
y = x*(sp.exp(1)-sp.exp(x))
y.subs({x:1/2})
def L(y):
return sp.simplify( sp.diff(y,x,x) - y )
f = L(y)
f

Out[2]:
$\displaystyle - e x - 2 e^{x}$
In [3]:
phi = x*(1-x)
a1 = sp.symbols('a1')
a1

Out[3]:
$\displaystyle a_{1}$
In [4]:
sp.expand(L(y)-f)

Out[4]:
$\displaystyle 0$
In [5]:
rho = sp.expand(L(a1*phi) - f)
rho

Out[5]:
$\displaystyle a_{1} x^{2} - a_{1} x - 2 a_{1} + e x + 2 e^{x}$
In [6]:
import numpy as np
yf = sp.lambdify(x,y,'numpy')
phif = sp.lambdify(x,phi,'numpy')
rhof = sp.lambdify((a1,x),rho,'numpy')
rhof(1.1,np.linspace(0,1,5))

Out[6]:
array([-0.2       ,  0.84137129,  2.18158346,  3.8664614 ,  5.95484549])
In [7]:
import matplotlib.pyplot as plt

In [8]:
xa = np.linspace(0,1,501)
plt.plot(xa,yf(xa),'k')
plt.plot(xa,phif(xa)*2.2);
#for a1n in np.linspace(-1,1,11):
#    plt.plot(xa,rhof(a1n,xa))

In [9]:
eqn = sp.integrate( rho*phi, (x,0,1) )
eqn

Out[9]:
$\displaystyle - \frac{11 a_{1}}{30} - \frac{23 e}{12} + 6$
In [10]:
a1star = sp.solve(eqn,a1)[0]
a1star

Out[10]:
$\displaystyle \frac{180}{11} - \frac{115 e}{22}$
In [11]:
galerkin_solution = a1star*phi
gsf = sp.lambdify(x,galerkin_solution,'numpy')
plt.plot(xa,yf(xa),'k')
plt.plot(xa,gsf(xa));

In [12]:
maxerrs = []
epsilons = np.linspace(-.01,.01,21)
for epsilon in epsilons:
maxerrs.append( np.abs((1+epsilon)*gsf(xa)-yf(xa)).max() )
plt.plot(epsilons,maxerrs,'o-')
plt.axvline([0])
plt.title('max error as coefficient varied')

Out[12]:
Text(0.5, 1.0, 'max error as coefficient varied')
In [16]:
plt.fill_between(xa,0*xa,rhof(float(a1star),xa)*phif(xa))
plt.grid()
plt.title('residual times phi');

In [ ]:


In [ ]: