from resources306 import *  # for expressionplot

Make up a BVP where we know the solution

x = sp.symbols('x')
# choose the solution in advance
y = sp.sin(sp.pi*x)*sp.exp(x)     # choose anything that satisfies y(0)=0, y(1)=0
expressionplot(y,x,0,1)

# Now choose a differential differential operator
def L(y):
    return sp.simplify(sp.diff(y,x,x) + 2*x*sp.diff(y,x) + 5*sp.exp(x)*y)

# and then figure our what f has to be in order for y to be the solution
f = L(y)
f

$\displaystyle \left(2 x \left(\sin{\left(\pi x \right)} + \pi \cos{\left(\pi x \right)}\right) + 5 e^{x} \sin{\left(\pi x \right)} - \pi^{2} \sin{\left(\pi x \right)} + \sin{\left(\pi x \right)} + 2 \pi \cos{\left(\pi x \right)}\right) e^{x}$

Now pretend we don't know the solution of L[y]=f and see if we can construct the optimal approximation using the Galerkin method

# basis of dimension 1
phi = x*(x-1)
a1_val = -6   # Pranav's suggestion
expressionplot( a1_val*phi, x, 0, 1,label=f'{a1_val} $\phi$')
expressionplot( y,x,0,1,label='solution y')
plt.legend();

We seek to optimize the approximation by choosing $a_1$ such that the residual is orthogonal to $\phi_1$

# compute the residual
a1 = sp.symbols('a1')
u = a1*phi
u

$\displaystyle a_{1} x \left(x - 1\right)$

rho = L(u)-f  # the residual
rho

$\displaystyle a_{1} \left(5 x \left(x - 1\right) e^{x} + 2 x \left(2 x - 1\right) + 2\right) - \left(2 x \left(\sin{\left(\pi x \right)} + \pi \cos{\left(\pi x \right)}\right) + 5 e^{x} \sin{\left(\pi x \right)} - \pi^{2} \sin{\left(\pi x \right)} + \sin{\left(\pi x \right)} + 2 \pi \cos{\left(\pi x \right)}\right) e^{x}$

# need rho to be orthogonal to each phi (of which there's only one here)
# obj = <rho,phi>
obj = sp.simplify(sp.integrate( rho*phi, (x,0,1)))
obj  # want this to be zero

$\displaystyle \frac{- 85665 \pi^{10} a_{1} - 5711 \pi^{12} a_{1} - 496857 \pi^{8} a_{1} - 1399195 \pi^{6} a_{1} - 1987428 \pi^{4} a_{1} - 1370640 \pi^{2} a_{1} - 365504 a_{1} + 134400 e a_{1} + 504000 e \pi^{2} a_{1} + 730800 e \pi^{4} a_{1} + 514500 e \pi^{6} a_{1} + 182700 e \pi^{8} a_{1} + 2100 e \pi^{12} a_{1} + 31500 e \pi^{10} a_{1} - 9480 e \pi^{7} - 840 e \pi^{9} - 1260 \pi^{9} - 11400 \pi^{7} - 33120 e \pi^{5} - 60 \pi^{11} - 19140 \pi^{5} - 24960 e \pi^{3} - 900 \pi^{5} e^{2} - 2700 \pi^{3} e^{2} - 1200 \pi e^{2} + 74640 \pi + 38400 e \pi + 26820 \pi^{3} + 1500 \pi^{7} e^{2} + 900 \pi^{9} e^{2}}{30 \left(64 + 240 \pi^{2} + 348 \pi^{4} + 245 \pi^{6} + 87 \pi^{8} + \pi^{12} + 15 \pi^{10}\right)}$

a1_optimal = sp.solve(obj,a1)[0]
display(a1_optimal)
float(a1_optimal)

$\displaystyle \frac{60 \pi \left(- 158 e \pi^{6} - 14 e \pi^{8} - 21 \pi^{8} - 190 \pi^{6} - 552 e \pi^{4} - \pi^{10} - 319 \pi^{4} - 416 e \pi^{2} - 15 \pi^{4} e^{2} - 45 \pi^{2} e^{2} - 20 e^{2} + 1244 + 640 e + 447 \pi^{2} + 25 \pi^{6} e^{2} + 15 \pi^{8} e^{2}\right)}{- 31500 e \pi^{10} - 2100 e \pi^{12} - 182700 e \pi^{8} - 514500 e \pi^{6} - 730800 e \pi^{4} - 504000 e \pi^{2} - 134400 e + 365504 + 1370640 \pi^{2} + 1987428 \pi^{4} + 1399195 \pi^{6} + 496857 \pi^{8} + 5711 \pi^{12} + 85665 \pi^{10}}$

$\displaystyle -4.54622036078795$

a1_val = float(a1_optimal)
expressionplot( a1_val*phi, x, 0, 1,label=f'{a1_val} $\phi$')
expressionplot( y,x,0,1,label='solution y')
plt.title('Galerkin approximation (blue)')
plt.legend();

Looks unimpressive, but ... let's look at the graph of $\rho \phi$:

expressionplot( rho.subs({a1:a1_optimal})*phi, x, 0,1)
plt.axhline(0,color='k',alpha=0.5)
plt.title('rho phi');

... it looks like the Galerkin orthogonality is satisfied

and in fact

sp.integrate( rho.subs({a1:a1_optimal})*phi, (x,0,1) )

$\displaystyle 0$

That the Galerkin approximation looks worse than Pranav's choice ($a_1=-6$) highlights the observation that the Galerkin approximation does not minimize $||e||^2 = \int_0^1 e^2$.