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:


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)$ → final values

  3. Find solution at final point $x_f$

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

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.

  • $\begingroup$ Is there a specific reason why you don't want to use a finite difference scheme? $\endgroup$
    – VoB
    Apr 24 '20 at 18:46
  • $\begingroup$ I'm unsure what you mean by this. I have only seen the shooting method used for Boundary Value questions. $\endgroup$
    – GavinK14
    Apr 24 '20 at 20:55
  • $\begingroup$ 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. $\endgroup$ Apr 25 '20 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 of y(x) for different values of y(-1)


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


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$")


  • 1
    $\begingroup$ Possible improvement for the OP: using a Newton method $\endgroup$
    – VoB
    Apr 28 '20 at 11:24
  • $\begingroup$ 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. $\endgroup$ Apr 28 '20 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));

    for k in range(0,N):

    return [X,Y,Y_prime];


while (b-a)>.0000001:
    if z>0:



Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.