请在下面找到针对您的问题的 Python 中Runge-Kutta 2 方法的实现。对于给定的（固定）和（变化）到的微分方程的积分。$-1$$5$$y'(-1)=-1$$y(-1)$

如您所见，使用此方法设置会得到而对于它会产生。例如，要找到的值，您需要对和之间值进行二分搜索（或任何其他寻根算法）。但是请注意，根的最终精度将受到积分方法精度的限制，因此寻求更精确的解决方案将需要增加步数。$y(-1) = 1$$y'(5) \approx 1.34$$y(-1) = 2$$y'(5) \approx -5.63$$y'(5) = 0$$y(-1)$$1$$2$

阴谋：

代码：

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()

这是使用 Forward Euler 方法进行积分和二分法找到的最小解决方案。我发现与我使用 Mathematica 时得到的值相同。$y(-1)$$y(-1)\approx 1.1926$

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)