Численное интегрирование обыкновенных дифференциальных уравнений

Решить краевую задачу для дифференциального уравнения второго порядка

с краевыми условиями первого рода: y(x₀) = y₀, y(x1) = y1 – на отрезке x ∈ [x₀, x₁].

Исходные данные:

Для наглядности дальнейших действий запишем общий вид дифференциального уравнения:

y'' + p(x)y' + q(x)y = f(x)

Теорию рекомендую почитать здесь (стр. 54 - 56)!

В процессе решения будет получена трехдиагональная СЛАУ, которую можно решить методом прогонки или методом Гаусса. В моем коде решение осуществляется методом Гаусса, т. к. он точнее.

Прежде чем запускать основной код, нужно сформировать код для метода Гаусса и проинициализировать его.

Вспомогательный код для решения СЛАУ методом Гаусса (подробнее в этой статье):

import numpy

def Metod_Gaussa(arr, brr):
    for k in range(arr.shape[0] - 1):
        #поиск строки с максимальным элементом
        max_elem = 0
        str = 0
        for i in range (k, arr.shape[0]):
            if abs(arr[i,k]) > abs(max_elem):
                max_elem = arr[i,k]
                str = i
        #меняем местами строки квадратной матрицы
        change = numpy.repeat(arr[k], 1)
        arr[k], arr[str] = arr[str], change
        #меняем местами элементы вектора-столбца
        change = numpy.repeat(brr[k], 1)
        brr[k], brr[str] = brr[str], change
        #делим полученную строку на max_elem
        arr[k] = arr[k] / max_elem
        brr[k] = brr[k] / max_elem
        #домножаем строку на коэффициенты и вычитаем ее из остальных строк
        for i in range (k + 1, arr.shape[0]):
            factor = arr[i,k]
            arr[i] = arr[i] - arr[k] * factor
            brr[i] = brr[i] - brr[k] * factor

    #находим неизвестные
    arg = [brr[brr.shape[0] - 1] / (arr[arr.shape[0] - 1, arr.shape[0] - 1])]
    for i in range(arr.shape[0] - 2, -1, -1):
        n = brr[i]
        for j in range(len(arg)):
            n = n - arg[j] * arr[i, arr.shape[0] - 1 - j]
        arg.append(n)

    #переворачиваем значения в списке
    otv = []
    for i in reversed(arg): otv.append(i)
    return otv

Этот код следует добавить в отдельный файл в ту же папку PyCharm, что и основной код:

Основной код:

import math
import numpy
import matplotlib.pyplot as plt
from Project1 import Gaussa

#задаем функции для дифференциального уравнения вида y" + py' + qy = f
def p(x): return -5
def q(x): return 6
def f(x): return math.exp((3 * x)) + x
x = [0, 3]
y = [2, 4]
n = 60 #число интервалов

def reshenie(p, q, f, x, y, n):
    h = (x[1] - x[0]) / n #шаг
    N = n - 1
    A = numpy.zeros([N, N])
    A[0, 0] = h ** 2 * q(x[0] + h) - 2
    A[0, 1] = 1 + 0.5 * h * p(x[0] + h)
    A[N - 1, N - 2] = 1 - 0.5 * h * p(x[1] - h)
    A[N - 1, N - 1] = h ** 2 * q(x[1] - h) - 2
    xi = x[0] + h
    for i in range(1, N - 1):
        xi += h
        A[i, i - 1] = 1 - 0.5 * h * p(xi)
        A[i, i] = h ** 2 * q(xi) - 2
        A[i, i + 1] = 1 + 0.5 * h * p(xi)
    B = numpy.zeros([N])
    B[0] = h ** 2 * f(x[0] + h) - y[0] * (1 - 0.5 * h * p(x[0] + h))
    B[N - 1] = h ** 2 * f(x[1] - h) - y[1] * (1 + 0.5 * h * p(x[1] - h))
    xi = x[0] + h
    global mas_xi
    mas_xi = [xi]
    for i in range(1, N - 1):
        xi += h
        mas_xi.append(xi)
        B[i] = h ** 2 * f(xi)
    mas_xi.append(x[1] - h)
    #с помощью метода Гаусса находим yi
    global yi
    yi = []
    for j in Gaussa.Metod_Gaussa(A, B):
        yi.append(float(j))
    point = [[] for i in range(len(yi))]
    for i in range(len(yi)):
        point[i].append(mas_xi[i])
        point[i].append(yi[i])
    return point

print(reshenie(p, q, f, x, y, n))

#график
mas_xi.insert(0, x[0]), yi.insert(0, y[0])
mas_xi.append(x[1]), yi.append(y[1])
plt.plot(mas_xi, yi)
plt.grid()
plt.show()

Результат:

Теперь разберемся с тем, что получили. Для этого рассмотрим пример попроще.

Далее применим тот же код, только в самом верху заменим параметры:

#задаем функции для дифференциального уравнения вида y" + py' + qy = f
def p(x): return x
def q(x): return -2
def f(x): return 2
x = [0, 4]
y = [0, 16]

Результат - парабола, что и следовало ожидать: