Solving a 2nd Order ODE using Finite Difference Method when Mixed Boundary Conditions are given

Question

The problem I'm looking at is

$$y'' + 3.05 y' -2.85 = 0 $$

with the boundary conditions $y(0) = 1$ and $y'(1) = 0.0305$.

After obtaining the algebraic set of equations using FDM, I'm not sure how the second condition would fit in to solve them? Please guide me on how we can best approach this problem.

P.S. We can solve the system in a straightforward way if $y(1)$ value was given instead of $y'(1)$. Is there a procedure to convert this Neumann conditon to Dirichlet condition (without the unknown integration constant)?

Answer 1

You could use a leapfrog scheme. Make the node number even, segment number odd, $x_0=a$, $x_{2N+1}=b$. Compute $y$ values at even nodes and $v$ values at odd nodes. $v$ is an approximation of $y'$. As finite difference approximation of the DE use \begin{align} y_{2k+2}-y_{2k}&=2hv_{2k+1}\\ v_{2k+1}-v_{2k-1}&=2h·2.85-h·3.05·(v_{2k+1}+v_{2k-1}) \end{align} Then the boundary conditions can be directly implemented as $y_0=1$ and $v_{2N+1}=0.0305$.

If the state vector is implemented "naturally" by implementing $y$ and $v$ values in the index order, this gives a system with tri-diagonal system matrix.

In code

a=3.05; r=2.85
y0=1; v1=0.0305

# --- exact solution ---
d = r/a; c = np.exp(a)*(v1-d)/a
y_ex=lambda x: y0+c*(1-np.exp(-a*x)) + d*x

# --- use sparseness of tridiagonal matrix ---
from scipy import sparse

N=15  # indices 0,...,2N+1, 2N+2 nodes, 2N+1 subintervals
x = np.linspace(0,1,2*N+2)
h = x[1]-x[0]

# equations from v_1 to y_{2N} as central variable
L = np.array([  1,a*h-1]*N)
D = np.array([2*h,   0 ]*N)
R = np.array([ -1,a*h+1]*N)
b = np.array([  0,2*h*r]*N)

A = sparse.diags([L,D,R],[0,1,2],shape=(2*N,2*N+2)).tocsr()

yv = np.zeros(2*N+2)
yv[0]=y0; yv[-1]=v1
b -= A.dot(yv)

yv[1:-1] = sparse.linalg.spsolve(A[:,1:-1],b)

# --- plot ---
plt.plot(x[::2], yv[::2],x,y_ex(x))
plt.grid(); plt.show()

This gives the plot

where the difference is barely recognizable. With smaller grid steps, like for $N=50$, no difference is visible.

Answer 2

Another possibility would be using (second order accurate) central finite differences for the inner points, and one-sided second order accurate FD at the boundary. Then the "Neumann-like" boundary condition $y'(1) = 0.0305$ translates into $$ \frac{1.5y_N -2 y_{N-1} +0.5 y_{N-2} }{h} = 0.0305$$ and can be integrated into the linear system that arises from applying finite differences for the other gridpoints.

clear;
close all;
clc;

N = 100;
h = 1/(N - 1);

% Use for both y'' and y' central finite differences
Diag_Elem = -2/h^2;
Lower_Elem = -3.05/(2*h) + 1/h^2;
Upper_Elem =  3.05/(2*h) + 1/h^2;

A = gallery('tridiag', Lower_Elem * ones(N-1, 1), ... 
                       Diag_Elem  * ones(N, 1), ...
                       Upper_Elem * ones(N-1, 1) );
                     
% Implement the y(0) = 1 B.C.:
A(1, 1) = 1;
A(1, 2) = 0;
                                   
% Implement the y'(1) = 0.0305 B.C.:
A(N, N-2) =  1/(2*h);
A(N, N-1) = -2/(h);
A(N, N)   =  3/(2*h);

% Fill RHS
b = 2.85 * ones(N, 1);
b(1) = 1;
b(N) = 0.0305;

y_numeric = A\b;

%% Analytical solution (see: https://www.wolframalpha.com/input?i=y%27%27+%2B+3.05+y%27+-2.85+%3D+0 )
c_1 = -19.08670023;
c_2 = 1 + 0.327869 * c_1;

x = linspace(0, 1, N);

y_analytic = -0.327869 * c_1 * exp(-3.05 * x) + c_2 + 0.934426 * x;

%% Plot 
plot(x, y_analytic, x, y_numeric, 'LineWidth', 2.0);
legend('Analytic', 'Numeric');
xlabel('$x$','Interpreter','latex', 'Fontsize', 18);
ylabel('$y$','Interpreter','latex', 'Fontsize', 18);
title('$\frac{d^2y}{dx^2} + 3.05 \frac{dy}{dx} -2.85 = 0$','Interpreter','latex', 'Fontsize', 22);