x = np.linspace(0, 1, 100);
xe = np.linspace(0, 1, 6);
ep = 0.1;
m1 = (-1 + (1+ep)**0.5)/ep; m2 = (-1 - (1+ep)**0.5)/ep; A = 1/(np.exp(m1) - np.exp(m2)); B = -A;
ye = A*np.exp(m1*xe) + B*np.exp(m2*xe);
plt.plot(xe, ye,'sr' ,label='Exact'); plt.xlabel('x'); plt.ylabel('y(x)');
yo0 = np.exp((x-1)/2); yo1 = 1/2*(1-x)*np.exp((x-1)/2);
yo = yo0 + ep*yo1;
plt.plot(x, yo, '--k', label='Outer sol');
yi = np.exp(-1/2)*(1-np.exp(-2*x/ep));
plt.plot(x, yi, '--r', label='Inner sol');
from scipy.integrate import solve_bvp
def fun(x,y):
ep = 0.1;
return [y[1], -2*y[1]/ep+y[0]/ep]
def bc(ya, yb):
return [ya[0]-0.0, yb[0]-1.0]
xi = np.linspace(0,1,10)
y_a = np.zeros((2,xi.size))
res_a = solve_bvp(fun, bc, xi, y_a)
xplot = np.linspace(0,1,100);
yplot = res_a.sol(xplot)[0];
plt.plot(xplot, yplot, '-.k', label='Using SciPy');
plt.legend();