Solving a Boundary Value Question $\frac{d^2y}{dx^2}=y\cos(x)+\frac{\sin(x)}{x^2+2}$ using Python

I'm looking to solve this boundary value question using the shooting method!

$$\frac{d^2y}{dx^2}=y\cos(x)+\frac{\sin(x)}{x^2+2}$$ given the initial values:

$$y'(x=-1)=-1\\y'(x=5)=0$$

I'm aware of the steps I should follow

1. Guess unknown initial values $$v_i$$

2. Solve ODE with these values: $$f (x|v_i)$$ → ﬁnal values

3. Find solution at ﬁnal point $$x_f$$

4. Solve $$f (x_f |v_i) − y_f = 0$$root ﬁnding!

I'm quite new to python 3.7 so if someone could help me code this problem or provide me with some tips/hints I would really appreciate it.

• Is there a specific reason why you don't want to use a finite difference scheme? – VoB Apr 24 at 18:46
• I'm unsure what you mean by this. I have only seen the shooting method used for Boundary Value questions. – GavinK14 Apr 24 at 20:55
• The shooting method easiest to implement would be probably also finite-difference based, so it is not the choice of finite-difference vs. shooting. The distinction is in using the shooting method vs. formulating the problem as a linear algebra problem and using a standard linear algebra package for it. – Maxim Umansky Apr 25 at 16:29

Please find below an implementation of the Runge-Kutta 2 method in Python for your problem. This takes care of integrating the differential equation from $$-1$$ to $$5$$ for a given value of $$y'(-1)=-1$$ (fixed) and $$y(-1)$$ (changing).

As you can see, setting $$y(-1) = 1$$ gives $$y'(5) \approx 1.34$$ with this method and for $$y(-1) = 2$$ it yields $$y'(5) \approx -5.63$$. To find the value for which $$y'(5) = 0$$ you would for instance need to do a dichotomic search on the value of $$y(-1)$$ between $$1$$ and $$2$$ (or any other root-finding algorithm). Note however that the ultimate precision on the root will be limited by the precision of the integration method, so going for a more precise solution would require to increase the number of steps.

Plot:

Code:


import numpy as np

xmin = -1
xmax = 5
Num_points = 600 #total number of steps
dx = (xmax-xmin)/Num_points
X = np.linspace(xmin, xmax, Num_points+1)
dy_min = -1 #value of dy/dx at xmin, given by the problem

def RK2_method(y_min):
y_list = np.zeros(Num_points+1)
dy_list = np.zeros(Num_points+1)
y_list[0] = y_min
dy_list[0] = dy_min
for k in range(Num_points):
y_half_step = y_list[k] + dx*dy_list[k]/2 #evaluating y and dy at n+1/2 according to the RK2 method
dy_half_step = dy_list[k] + dx/2*(y_list[k]*np.cos(X[k]) + np.sin(X[k])/(X[k]**2+2))
new_y = y_list[k] + dx*dy_half_step
new_dy = dy_list[k] + dx*(y_half_step*np.cos(X[k]+dx/2) + np.sin(X[k]+dx/2)/((X[k]+dx/2)**2+2))
y_list[k+1] = new_y
dy_list[k+1] = new_dy
return y_list, dy_list

y_list_1, dy_list_1 = RK2_method(y_min=1)
y_list_2, dy_list_2 = RK2_method(y_min=2)

print(dy_list_1[-1]) #dy/dx(5) for y(-1) = 1 is > 0
print(dy_list_2[-1]) #dy/dx(5) for y(-1) = 2 is < 0

#####Plotting#####

import matplotlib.pyplot as plt
plt.rc('font', size=24)

fig, ax = plt.subplots(1)
fig.suptitle(r"Solving $$\frac{d^2y}{dx^2} = y \cos(x) + \frac{\mathrm{\sin(x)}}{x^2+2}, y'(-1)=-1$$")

ax.plot(X, y_list_2, 'r-', lw=2, label=r"$$y(-1)=2$$")
ax.plot(X, y_list_1, 'b-', lw=2, label=r"$$y(-1)=1$$")

ax.set_xlabel(r'$$x$$')
ax.set_ylabel(r'$$y(x)$$')
ax.set_xlim(-1,5)
plt.legend(loc='best')
plt.show()


• Possible improvement for the OP: using a Newton method – VoB Apr 28 at 11:24
• That is correct, as the Newton method would allow to find the root faster. At the same time, the precision of the root finding is ultimately limited by the precision of the integration method, so depending on the target precision, it might not be that much longer to go for a dichotomic search instead. – QuantumApple Apr 28 at 11:33

Here's a minimal solution using Forward Euler method for integration and dichotomy to find $$y(-1)$$. I found that $$y(-1)\approx 1.1926$$ which is the same value that I got when I used Mathematica.

from math import *;
import matplotlib.pyplot as plt

def solve(a):
#Solve the IVP y''(x)=f(y(x),x) with the initial conditions y(-1)=a and y'(-1)=1 using Forward Euler method.
N=30000;dx=6/N;X=list(-1+k*dx for k in range(0,N+1));
Y=list(0 for k in range(0,N+1));
Y_prime=list(0 for k in range(0,N+1));
Y[0]=a;Y_prime[0]=-1;

for k in range(0,N):
x=X[k];
up=Y_prime[k];
vp=Y[k]*cos(x)+sin(x)/(x*x+2);
Y[k+1]=Y[k]+dx*up;
Y_prime[k+1]=Y_prime[k]+dx*vp;

return [X,Y,Y_prime];

a=1;b=2;

while (b-a)>.0000001:
c=(a+b)/2;
[X,Y,Y_prime]=solve(c);
z=Y_prime[-1];
if z>0:
a=c;
else:
b=c;

print(a)

plt.plot(X,Y)
plt.plot(X,Y_prime)