5. Инженерные расчеты python

1. Научные и Инженерные расчеты

НАУЧНЫЕ И
ИНЖЕНЕРНЫЕ
РАСЧЕТЫ

2. Библиотека SciPy

SciPy — это библиотека Python с открытым исходным кодом,
предназначенная для решения научных, инженерных и
математических проблем.
Она построена на базе NumPy и позволяет управлять данными, а
также визуализировать их с помощью разных высокоуровневых
команд. Если вы импортируете SciPy, то NumPy отдельно
импортировать не нужно.

3. Библиотека SciPy

В SciPy есть набор пакетов для разных научных вычислений:
• cluster – алгоритмы кластерного анализа
• constants – физические и математические константы
• fftpack – быстрое преобразование Фурье
–
численное интегрирование,
дифференциальных уравнений и пр.
• integrate
• interpolate – интерполяция и сглаживание сплайнов
решение

4. Библиотека SciPy

• io – ввод и вывод
• linalg – линейная алгебра
• ndimage – N-размерная обработка изображений
• odr – метод ортогональных расстояний
• optimize – оптимизация и численное решение уравнений
• signal – обработка сигналов
• sparse – разреженные матрицы

5. Библиотека SciPy

• spatial – разреженные структуры данных и алгоритмы
• special – специальные функции
• stats – статистические распределения и функции

6. SciPy.constants

from scipy.constants import *
print('g', g)
print('electron mass', electron_mass)
print(physical_constants['electron volt'])
print('nano', nano)

7. SciPy.constants

Углы

8. SciPy.constants

Время
Длина

9. SciPy.constants

Давление
Скорость

10. SciPy.constants

Сила, энергия,
мощность

11. Пример использования констант

Построить график кинетической энергии частицы, попадающей в
защитный экран космического корабля в зависимости от скорости
столкновения в долях скорости света. Энергию частицы выразить
в тротиловом эквиваленте.
Кинетическая энергия частицы сферической
формы радиуса R с плотностью материала ρ,
движущейся со скоростью

12. Пример использования констант

import numpy as np
import scipy.constants as CONST
import matplotlib.pyplot as plt
R = 0.0005
rho = 1000.0
m = rho * (4.0 / 3.0) * np.pi * R ** 3
# Массив долей скорости света
k = np.linspace(0 , 0.1 , 9)
print("k:", k)
# Энергия (массив)
Energy = 0.5 * m * (k * CONST.c) ** 2
# Энергия в тоннах ТНТ
E_ton_TNT = Energy / CONST.ton_TNT
# Энергия в килограммах ТНТ
E_kg_TNT = 1000.0 * E_ton_TNT
print("E_kg_TNT:", E_kg_TNT)
plt.figure(figsize = (6 , 3.5))
plt.plot (k , E_kg_TNT , 'r')
plt.xlabel('Скорость частицы, v/c')
plt.ylabel('Энергия частицы , кг ТНТ')
plt.legend(['m = {:3.2f} мг '.
format(m*1e6)])
plt.grid()
plt.show()

13. SciPy.interpolate

Одномерная интерполяция — это область построения
кривой, которая бы полностью соответствовала набору
двумерных точек данных. В SciPy есть функция interp1d,
которая
используется
для
создания
одномерной
интерполяции.

14. Интерполяция и аппроксимация

Дано:
• Набор экспериментальных точек (xi, yi), i = 1 .. n
• Часто точки зашумлены или получены с погрешностью
Цель:
• Построить функцию f(x), которая:
• Хорошо описывает зависимость между x и y
• Позволяет предсказывать значения в новых точках
• Учитывает природу данных
Два подхода:
• Интерполяция
• Аппроксимация

15. Интерполяция

Интерполяция – метод построения функции f(x), которая точно проходит
через все заданные точки (xi,yi):
Где применяется:
• Восстановление пропущенных значений в таблицах
• Построение гладких кривых по ключевым точкам (CAD, анимация)
• Приближение сложных функций при известных точных значениях
Ограничения:
• Чувствительна к шуму (будет проходить через выбросы)
• При большом количестве точек → сильные колебания

16. Основные методы интерполяции

Метод
Линейная
Формула
Особенности
отрезки прямых между Простая, негладкая
точками
(разрывы производной)
Полиномиальная
(Лагранж)
полином n−1 степени
Гладкая, но сильные
колебания при n>10
Сплайновая (Bсплайны, кубические)
кусочные полиномы
Гладкая, устойчивая,
оптимальная
учёт производных
Ещё более гладкая,
требует доп. данных
Эрмитова

17. SciPy.interpolate

interp1d (x, y, kind='linear' , axis =-1 , copy = True ,
bounds_error = None, fill_value = nan , accept_sorted = False)
kind – строка или целое число, указывающие тип интерполяции. Если
параметр задается в виде строки, то он может принимать одно из значений:
'linear', 'nearest', 'nearest-up', 'zero', 'slinear', 'quadratic', 'cubic', 'previous' или
'next'. Значения 'zero', 'slinear', 'quadratic' и 'cubic' относятся к сплайнинтерполяции нулевого, первого, второго или третьего порядка. Если
значения параметра равны 'previous' или 'next', то возвращается
предыдущее или следующее значение функции y; 'nearest-up' и 'nearest'
различаются при интерполяции нецелых чисел (например, 0,5, 1,5) тем, что
'nearest-up' округляет значения с избытком, а 'nearest' – с недостатком.

18. SciPy.interpolate

interp1d (x, y, kind='linear' , axis =-1 , copy = True ,
bounds_error = None, fill_value = nan , accept_sorted = False)
axis – целочисленная переменная. Определяет ось, вдоль которой
выполняется интерполяция. По умолчанию используется последняя ось;
copy – логическая переменная. Если значение параметра равно True,
создаются внутренние копии x и y. Если значение параметра равно False,
используются ссылки на x и y;

19. SciPy.interpolate

interp1d (x, y, kind='linear' , axis =-1 , copy = True ,
bounds_error = None, fill_value = nan , accept_sorted = False)
bounds_error – булева переменная. При значении переменной,
равном True, возникает ошибка ValueError всякий раз, когда
предпринимается попытка вычисления значения функции в точке,
выходящей за пределы диапазона изменения переменной x. Если значение
переменной bounds_error равно False (или 0), результатом работы функции
является NaN, если только значение параметра fill_value не равно
"extrapolate". Если значение параметра не задано, то возникает ошибка,
если только значение fill_value не равно "extrapolate";

20. SciPy.interpolate

interp1d (x, y, kind='linear', axis =-1 , copy = True ,
bounds_error = None, fill_value = nan, accept_sorted = False)
fill_value – параметр, использующийся при решении задачи экстраполяции.
Если fill_value="extrapolate", то решение будет получено и в случае, когда
значение аргумента, при котором вычисляется значение функции,
находится за пределами интервала [min(x), max(x)];
accept_sorted – логическая переменная. Если accept_sorted=False, значения
x могут быть расположены в любом порядке, и затем они автоматически
сортируются. Если значение параметра равно True, x должен быть массивом
монотонно возрастающих значений.

21. Пример интерполяции

Функция задана в виде таблицы
Вычислить значение функции в точке х* = 4,5. Решить задачу
интерполяции с помощью функции interp1d(), используя
линейную, квадратичную и кубическую интерполяцию. Дать
графическое решение задачи.

22. Пример интерполяции

from scipy.interpolate import interp1d
import numpy as np
import matplotlib.pyplot as plt
x = np.array([1, 2, 3, 4, 5])
y = np.array([0.1, 2.8, 2.3, 4.5, 5.3])
px = 4.5
f1 = interp1d(x, y, kind='linear')
xnew = np.linspace(1, 5, 101)
ynew1 = f1(xnew)
py1 = f1(px)
f2 = interp1d(x, y, kind='quadratic')
ynew2 = f2(xnew)
py2 = f2(px)
f3 = interp1d(x, y, kind='cubic')
ynew3 = f3(xnew)
py3 = f3(px)
plt.plot(xnew, ynew1, c='salmon', label='linear')
plt.plot(xnew, ynew2, c='steelblue', label='quadratic')
plt.plot(xnew, ynew3, c='limegreen', label='cubic')
plt.plot(px, py1,'X', c='crimson', ms=10,
label=f'flinear({px})={py1:.2f}')
plt.plot(px, py2, 'o', c='b', ms=10,
label=f'fquadratic({px})={py2:.2f}')
plt.plot(px, py3, 'X', c='g', ms=10,
label=f'fcubic({px})={py3:.2f}')
plt.plot(x, y, 'D', c='darkorchid', ms=8, label='data')
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()

23. SciPy.interpolate

В последнее время в вычислительной практике широкое
распространение получили B-сплайны (Basic spline). Они
используются как для интерполяции функций, так и в
качестве базисных функций при построении методов типа
конечных элементов. Метод конечных элементов – это
численный метод решения дифференциальных уравнений с
частными производными, а также интегральных уравнений,
возникающих, в частности, при решении задач прикладной
физики.

24. SciPy.interpolate

Для интерполяции с помощью B-сплайнов может быть
использована функция splrep(). В отличие от обычной
сплайн-интерполяции, «сшивка» элементарных B-сплайнов
производится не в точках (xi, yi), а в других точках,
координаты которых в функции splrep() могут либо
вычисляться
автоматически,
либо
задаваться
пользователем.

25. SciPy.interpolate

splrep(x, y, w=None, xb=None, xe=None, k=3, task=0,
s=None, t=None, full_output=0, per=0, quiet=1).
w – строго положительный массив весов той же длины, что и x с y.
Веса
используются
при
вычислении
коэффициентов
«взвешенного» сплайна методом наименьших квадратов. Если
погрешности в значениях y имеют стандартное отклонение,
заданное вектором d, то вес w должен быть равным 1/d. Значение
по умолчанию ones(len(x)): все весаравны 1;

26. SciPy.interpolate

splrep(x, y, w=None, xb=None, xe=None, k=3, task=0, s=None, t=None,
full_output=0, per=0, quiet=1).
xb, xe – переменная типа float. Интервал интерполяции. Если параметр опущен, по
умолчанию используются значения x[0] и x[-1] соответственно;
k – переменная типа int. Степень сплайна;
task – параметр, принимающий одно из значений {1, 0, -1}. Если task=0, находятся t и c
для заданного коэффициента сглаживания s. Если task=1, вычисляются t и c для
предыдущего значения коэффициента сглаживания. До этого должен быть осуществлен
вызов функции splrep() с task=0 или task=1 для того же набора данных. Если task=-1,
находится «взвешенный» сплайн для заданного набора узлов t;

27. SciPy.interpolate

splrep(x, y, w=None, xb=None, xe=None, k=3, task=0,
s=None, t=None, full_output=0, per=0, quiet=1).
s – переменная типа float. Коэффициент сглаживания. Величина
гладкости обеспечивает выполнение условий sum((w*(yg))**2,axis=0)≤ s, где g(x) – результат интерполяция набора
данных (x,y). По умолчанию s=m-sqrt(2*m), если указаны веса.
Здесь m – число точек таблицы. Значение параметра s=0, если вес
не указан;
t – массив, с помощью которого задаются узлы в случае, когда
значение параметра task=-1;

28. SciPy.interpolate

splrep(x, y, w=None, xb=None, xe=None, k=3, task=0,
s=None, t=None, full_output=0, per=0, quiet=1).
full_output – булева переменная. Если значение отлично от 0,
выводятся дополнительные выходные параметры функции
splrep();
per – булева переменная. Если значение не 0, исходная функция
считается периодической с периодом x[m-1] - x[0];
quiet – булева переменная. Если значение не равно 0,
подавляется вывод необязательных сообщений.

29. SciPy.interpolate

После построения сплайна вычислить его значения в заданных
точках можно с помощью функции splev(). Ее синтаксис:
splev(x, tck, der=0, ext=0)
x – массив точек, в которых нужно вычислить значение
сглаженного сплайна или его производных;
tck – кортеж или объект B-Spline. Если tck – кортеж, то это
последовательность длиной 3, возвращаемая функциями splrep()
или splprep(). tck содержит узлы, коэффициенты и степень
сплайна;

30. SciPy.interpolate

splev(x, tck, der=0, ext=0)
der – значение типа int. Порядок производной сплайна для
вычисления ее значений в заданных точках х (должен быть
меньше или равен степени сплайна k);
ext – значение типа int, задающее порядок действий в случае,
когда значение х выходит за пределы [x[0], x[-1]]. В зависимости
от
значения
параметра
ext
функция
возвращает
экстраполированное значение (ext=0), 0 (ext=1), сообщение
ValueError() (ext=2), граничное значение (y[0] или y[-1]) (ext=3).
Значение параметра ext по умолчанию – 0.

31. Пример интерполяции

По данным таблицы построить графики B-сплайнов 2-й и 3-й
степени на отрезке [1, 5]:

32. Пример интерполяции

import numpy as np
import matplotlib.pyplot as plt
from scipy.interpolate import splrep,
splev
x = np.array([1, 2, 4, 5])
y = np.array([8, 12, 6, 3])
x_smooth = np.linspace(1, 5, 200)
tck_2 = splrep(x, y, k=2)
y_smooth_2 = splev(x_smooth, tck_2)
tck_3 = splrep(x, y, k=3)
y_smooth_3 = splev(x_smooth, tck_3)
plt.plot(x_smooth, y_smooth_2, 'b-',
label='B-сплайн 2-й степени (k=2)')
plt.plot(x_smooth, y_smooth_3, 'g-',
label='B-сплайн 3-й степени (k=3)')
plt.plot(x, y, 'ro', label='Исходные
точки')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()

33. Аппроксимация

Аппроксимация
–
метод
построения
функции
которая минимизирует суммарную ошибку между точками и кривой
Где применяется:
• Обработка экспериментальных данных с шумом
• Машинное обучение (регрессия)
• Физические измерения (всегда есть погрешность)
Преимущества:
• Устойчива к выбросам
• Даёт гладкую, "осмысленную" кривую
• Позволяет управлять сложностью модели
f(x),

34. Оценка качества моделей

35. Оценка качества моделей

36. SciPy.optimize

Пакет scipy.optimize включает в себя реализацию следующих функций:
1.Условной и безусловной минимизации скалярных функций нескольких
переменных (minim) с помощью различных алгоритмов (симплекс НелдераМида, BFGS, сопряженных градиентов Ньютона, COBYLA и SLSQP)
2.Глобальной оптимизации (например: basinhopping, diff_evolution)
3.Минимизация остатков МНК (least_squares) и алгоритмы подгонки кривых
нелинейным МНК (curve_fit)
4.Минимизации скалярной функций одной переменной (minim_scalar) и
поиска корней (root_scalar)
5.Многомерные решатели системы уравнений (root) с использованием
различных алгоритмов (гибридный Пауэлла, Левенберг-Марквардт или
крупномасштабные методы, такие как Ньютона-Крылова).

37. SciPy.optimize

Функция scipy.optimize.curve_fit используется для аппроксимации
данных с помощью заданной функции. Она подбирает параметры
функции таким образом, чтобы минимизировать сумму квадратов
разностей между предсказанными значениями функции и исходными
данными.

38. SciPy.optimize

scipy.optimize.curve_fit(func, xdata, ydata, p0=None)
func: функция, которая будет использоваться для аппроксимации
данных. Должна принимать в качестве первого аргумента массив
xdata и остальные параметры, которые нужно подобрать.
Возвращает значения функции.
xdata: массив значений для независимой переменной.
ydata: массив значений для зависимой переменной.
p0: начальное приближение для параметров функции. Если не
указано, то используются значения по умолчанию.

39. SciPy.optimize

scipy.optimize.curve_fit(func, xdata, ydata, p0=None)
Возвращает кортеж из двух элементов:
- popt: оптимальные значения параметров функции, которые
минимизируют ошибку аппроксимации.
- pcov: ковариационная матрица параметров.

40. Пример аппроксимации

Выполнить аппроксимацию табличных данных тремя заданными
моделями:
y ax bx c
2
y ae
bx
a
y
c
x b
c
Отсортировать
полученные
модели
детерминации (R2). Построить график
x
y
1
2,8
2
1,8
3
1,5
по
5
1,3
коэффициенту
7
1,1
8
1

41. Пример аппроксимации

import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit
x = np.array([1, 2, 3, 5, 7, 8])
y = np.array([2.8, 1.8, 1.5, 1.3, 1.1, 1])
def quadratic_model(x, a, b, c):
return a * x ** 2 + b * x + c
def exp_model(x, a, b, c):
return a * np.exp(-b * x) + c
def hyperbole_model(x, a, b, c):
return a / (b * x + c)
# Функция для вычисления R²
def calculate_r2(y_true, y_pred):
ss_res = np.sum((y_true - y_pred) ** 2)
ss_tot = np.sum((y_true - np.mean(y_true)) ** 2)
r2 = 1 - (ss_res / ss_tot)
return r2
models = {
'Квадратичная': quadratic_model,
'Экспоненциальная': exp_model,
'Гиперболическая': hyperbole_model
}
eq = {
'Квадратичная': '{:.3f}x²+{:.3f}x+{:.3f}',
'Экспоненциальная': '{:.3f}exp({:.3f}x)+{:.3f}',
'Гиперболическая': '{:.3f}/({:.3f}x+{:.3f})'
}

42. Пример аппроксимации

results = {}
for name, model in models.items():
popt, _ = curve_fit(model, x, y, maxfev=5000)
y_pred = model(x, *popt)
r2 = calculate_r2(y, y_pred)
results[name] = {
'params': popt,
'r2': r2,
'y_pred': y_pred,
'label': eq[name].format(*popt)
}

43. Пример аппроксимации

# Сравнительная таблица
print("="*54)
print("СРАВНЕНИЕ МОДЕЛЕЙ")
print("="*54)
print(f"|{'Модель':<16}|{'R²':<6}|{'Ранг'}|{'Уравнение':<23}|")
print("-"*54)
sorted_models = sorted(results.items(), key=lambda x: x[1]['r2'], reverse=True)
for rank, (name, result) in enumerate(sorted_models, 1):
print(f"|{name:<16}|{result['r2']:.4f}|{rank}
|{result['label']:<23}|")

44. Пример аппроксимации

# График
plt.scatter(x, y, label='Данные', c='k')
xnew = np.linspace(x.min(), x.max(), 200)
for name, result in results.items():
ynew = models[name](xnew, *result['params'])
plt.plot(xnew, ynew, label=result['label']+f'(R²={result["r2"]:.3f})')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.tight_layout()
plt.show()

45. Численное интегрирование

В инженерных расчетах довольно часто возникает необходимость
вычисления определенного интеграла от некоторой функции
b
f x dx
a
где f(x) – некоторая заданная на отрезке [a, b] непрерывная функция.

46. Численное интегрирование

Геометрический смысл интегрирования –
это нахождение площади фигуры,
ограниченной графиком функции f(x),
отрезком оси абсцисс, прямой x=a и
прямой x=b, т.е. вычисление интеграла
равносильно
вычислению
площади
криволинейной трапеции.

47. Численное интегрирование

Численное интегрирование (историческое название: (численная)
квадратура) – вычисление значения определённого интеграла (как
правило, приближённое) от функций, заданных либо в явном виде,
либо в виде таблицы.
К численному интегрированию прибегают тогда, когда точными
методами интеграл вычислить либо невозможно, либо сложно.

48. Численное интегрирование

Разобьем промежуток интегрирования на n участков шириной
b a
h
n
и заменим искомый интеграл на сумму интегралов (площадей фигур
под графиком функции) на полученных элементарных участках:
b
n
k
a
k 1 xk 1
f x dx f x dx
Заменяя подынтегральную функцию полиномами получаем формулы
для численного решения следующими методами:

49. Методы прямоугольников

На каждом участке интегрирования функция f(x) заменяется
полиномом нулевой степени – отрезком, параллельным оси абсцисс.
Различают метод левых, правых и центральных прямоугольников.

50. Метод левых прямоугольников

b
n 1
a
k 0
f x dx h f a kh

51. Метод правых прямоугольников

b
n
a
k 1
f x dx h f a kh

52. Метод центральных прямоугольников

h
f a kh
a f x dx h
2
k 1
b
n

53. Методы трапеций

В этом методе функция f(x) на каждом участке заменяется полиномом
первой степени.
f a f b n 1
f a kh
a f x dx h
2
k 1
b

54. Метод Симпсона

Подынтегральная функция f (x)
заменяется на каждом участке
полиномом второй степени –
параболой, проходящей через три
узла
n 1
n 2
h
f a kh 2 f a kh
a f x dx 3 f a f b 4k
1,3,5...
k 2,4,6...
b

55. SciPy.integrate

С помощью scipy.integrate можно вычислять определенные и
неопределенные интегралы одномерных и многомерных функций.
Также можно решать обыкновенные дифференциальные уравнения и
системы дифференциальных уравнений.
Для численного интегрирования одномерных функций можно
использовать функции quad, trapz, simpson, rromb. Для численного
интегрирования многомерных функций доступны функции dblquad,
tplquad, nquad.

56. SciPy.integrate

Методы для интегрирования функций с использованием объекта функции.
quad – Интегрирование общего назначения.
dblquad – Интегрирование общего назначения для двойного интеграла.
tplquad – Интегрирование общего назначения для тройного интеграла
fixed_quad – Интегрирование функции func(x) с использованием метода
Гауссовой квадратуры порядка n.

57. SciPy.integrate

Численные методы
интервалами:
интегрирования
с
фиксированными
trapezoid – по правилу трапеций.
cumulative_trapezoid – куммулятивное вычисление интеграла
по правилу трапеций.
simpson – интегрирование по методу Симпсона.
romb – Использование интегрирования Ромберга для
вычисления интеграла от (2^k) равномерно распределенных
отрезков.

58. Пример

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from scipy.interpolate import lagrange
from scipy import integrate
def left_rectangle(f, a, b, n):
h = (b - a) / n
x = np.linspace(a, b - h, n)
return h * np.sum(f(x))
def right_rectangle(f, a, b, n):
h = (b - a) / n
x = np.linspace(a + h, b, n)
return h * np.sum(f(x))
def center_rectangle(f, a, b, n):
h = (b - a) / n
x = np.linspace(a + h/2, b - h/2, n)
return h * np.sum(f(x))

59. Пример

def trapezoidal(f, a, b, n):
h = (b - a) / n
x = np.linspace(a, b, n+1)
y = f(x)
return h * (np.sum(y) - (y[0] + y[-1])/2)
def simpson(f, a, b, n):
# Если количество интевалов нечетное - увеличиваем на 1
if n % 2 != 0:
n += 1
h = (b - a) / n
x = np.linspace(a, b, n+1)
y = f(x)
return (h/3) * ((y[0] + y[-1]) + 4 * np.sum(y[1:-1:2]) + 2 * np.sum(y[2:-2:2]))

60. Пример

def draw_analytic_int(f, a, b):
x = np.linspace(a-0.05, b+0.05, 200)
y = f(x)
plt.axvline(x=0, c='gray', ls='-', lw=1)
plt.axhline(y=0, c='gray', ls='-', lw=1)
plt.plot(x, y, c='steelblue', lw=3)
plt.fill_between(x, y, where=(x>a) & (x<b), color='steelblue', alpha=0.3)
plt.grid(ls=':')
def print_result(i_num, i_an, name_method):
error = abs(i_num - i_an)
rel_error = error / abs(i_an)
print(f"| {name_method:<26} | {i_num:^9.6f} | {error:^8.6f} |
{rel_error:^9.2%}|")

61. Пример

a, b = 0, 2
n = 8
x = sp.symbols('x')
expr = x ** 2 + 0.5 * sp.sin(5 * x) + 0.2
i_analytic = float(sp.integrate(expr, (x, a, b)))
f = sp.lambdify(x, expr, "numpy")
print("=" * 64)
print(f"ЧИСЛЕННОЕ ИНТЕГРИРОВАНИЕ")
print(f"Функция: f(x) = x² + 0.5·sin(5x) + 0.2")
print(f"Интервал: [{a}, {b}]")
print(f"Количество отрезков: n = {n}")
print(f"Точное значение интеграла: {i_analytic:.6f}")
print("=" * 64)
print("РЕЗУЛЬТАТЫ РУЧНОЙ РЕАЛИЗАЦИИ:")
print("-" * 64)
print(f"| {'Метод':^26} | {'Значение':^9} | {'Ошибка':^8} |{'Отн.ошибка':^10}|")
print("-" * 64)

62. Пример

# Визуализация аналитического решения
plt.subplot(2, 3, 1)
draw_analytic_int(f, a, b)
plt.title('Аналитическое решение')
# Метод левых прямоугольников
i_left_rec = left_rectangle(f, a, b, n)
print_result(i_left_rec, i_analytic, 'Левые прямоугольники')
# Визуализация метода левых прямоугольников
x = np.linspace(a, b, n+1)
y = f(x)
plt.subplot(2, 3, 2)
draw_analytic_int(f, a, b)
for i in range(n):
X = [x[i], x[i], x[i+1], x[i+1]]
Y = [0, y[i], y[i], 0]
plt.fill(X, Y, c='gold', ec='peru', lw=2)
plt.title('Метод левых прямоугольников')

63. Пример

# Метод правых прямоугольников
i_right_rec = right_rectangle(f, a, b, n)
print_result(i_right_rec, i_analytic, 'Правые прямоугольники')
# Визуализация метода правых прямоугольников
plt.subplot(2, 3, 3)
draw_analytic_int(f, a, b)
for i in range(n):
X = [x[i], x[i], x[i+1], x[i+1]]
Y = [0, y[i+1], y[i+1], 0]
plt.fill(X, Y, c='gold', ec='peru', lw=2)
plt.title('Метод правых прямоугольников')

64. Пример

# Метод центральных прямоугольников
i_center_rec = center_rectangle(f, a, b, n)
print_result(i_center_rec, i_analytic, 'Центральные прямоугольники')
h=(b-a)/n
x_center = np.linspace(a, b, 2*n+1)
y_center = f(x_center)
# Визуализация метода центральных прямоугольников
plt.subplot(2, 3, 4)
draw_analytic_int(f, a, b)
for i in range(n):
X = [x_center[2*i], x_center[2*i], x_center[2*i+2], x_center[2*i+2]]
Y = [0, y_center[2*i+1], y_center[2*i+1], 0]
plt.fill(X, Y, c='gold', ec='peru', lw=2)
plt.title('Метод центральных прямоугольников')

65. Пример

# Метод трапеций
i_trapezoidal = trapezoidal(f, a, b, n)
print_result(i_trapezoidal, i_analytic,
'Трапеций')
# Визуализация метода трапеций
plt.subplot(2, 3, 5)
draw_analytic_int(f, a, b)
for i in range(n):
X = [x[i], x[i], x[i+1], x[i+1]]
Y = [0, y[i], y[i+1], 0]
plt.fill(X, Y, c='gold', ec='peru',
lw=2)
plt.title('Метод трапеций')
plt.tight_layout()
# Зум
n1, n2 = 2, 4
ax_zoom = plt.axes([0.41, 0.28, 0.12, 0.14])
draw_analytic_int(f, x[n1-1], x[n2+1])
for i in range(n1-1,n2+1):
X = [x[i], x[i], x[i+1], x[i+1]]
Y = [0, y[i], y[i+1], 0]
ax_zoom.fill(X, Y, c='gold', ec='peru',
lw=2)
ax_zoom.axis([x[n1]-0.05, x[n2]+0.05,
min(y[n1:n2+1])-0.1, max(y[n1:n2+1])+0.1])

66. Пример

# Метод Симпсона
i_simpson = simpson(f, a, b, n)
print_result(i_simpson, i_analytic, 'Симпсона')
n_simpson = n if n % 2 == 0 else n+1
x_simpson = np.linspace(a, b, n_simpson+1)
y_simpson = f(x_simpson)
# Визуализация метода Симпсона
plt.subplot(2, 3, 6)
draw_analytic_int(f, a, b)
for i in range(0, n_simpson, 2):
X_par = x_simpson[i:i+3]
Y_par = y_simpson[i:i+3]
polynomial = lagrange(X_par, Y_par)
x_plot = np.linspace(X_par[0], X_par[2], 100)
y_plot = polynomial(x_plot)
plt.fill_between(x_plot, y_plot, color='gold', ec='peru', lw=2)
plt.title('Метод Симпсона')

67. Пример

# Зум
n1, n2 = 1, 2
ax_zoom2 = plt.axes([0.74, 0.28, 0.12, 0.14])
draw_analytic_int(f, x_simpson[2*(n1-1)], x_simpson[2*(n2+1)])
for i in range(2*(n1-1), 2*(n2+1), 2):
X_par = x_simpson[i:i+3]
Y_par = y_simpson[i:i+3]
polynomial = lagrange(X_par, Y_par)
x_plot = np.linspace(X_par[0], X_par[2], 100)
y_plot = polynomial(x_plot)
ax_zoom2.fill_between(x_plot, y_plot, color='gold', ec='peru', lw=2)
ax_zoom2.axis([x_simpson[2*n1]-0.05, x_simpson[2*n2]+0.05,
min(y_simpson[2*n1:2*n2+1])-0.1, max(y_simpson[2*n1:2*n2+1])+0.1])
plt.show()

68. Пример

# ============== МЕТОДЫ SCIPY ==============
print("\n" + "=" * 64)
print("ФУНКЦИИ SCIPY.INTEGRATE:")
print("=" * 64)
result_quad, error_quad = integrate.quad(f, a, b)
rel_error_quad = abs(result_quad - i_analytic)/abs(i_analytic)
print(f"| {'quad':<26} | {result_quad:^9.6f} | {error_quad:^8.6f} |
{rel_error_quad:^9.2%}|")
n_gauss = 5 # количество узлов
result_fixed_quad = integrate.fixed_quad(f, a, b, n=n_gauss)[0]
print_result(result_fixed_quad, i_analytic,'fixed_quad')

69. Пример

# Методы scipy для табличных данных
x_scipy = np.linspace(a, b, n+1)
y_scipy = f(x_scipy)
result_trapezoid = integrate.trapezoid(y_scipy, x_scipy)
print_result(result_trapezoid, i_analytic,'trapezoid')
result_simpson = integrate.simpson(y_scipy, x_scipy)
print_result(result_simpson, i_analytic,'simpson’)
# Для romb количество интервалов должно быть степенью 2
result_romberg = integrate.romb(y_scipy, dx = (b-a)/n)
print_result(result_romberg, i_analytic,'romberg')
print("-" * 64)

70. Пример

71. Пример

from scipy import integrate
import numpy as np
import matplotlib.pyplot as plt
# Функция: f(x) = x² + sin(x)
def f(x):
return x**2 + np.sin(x)
# Интеграл от 0 до 2
result, error = integrate.quad(f, 0, 2)
print("="*50)
print("ПРИМЕР 1: scipy.integrate.quad")
print("="*50)
print(f"Функция: f(x) = x² + sin(x)")
print(f"Интеграл ∫₀² f(x) dx = {result:.8f}")
print(f"Оценка ошибки: {error:.2e}")

72. Пример

# Функция: f(x, y) = x * y
def f(y, x): # порядок аргументов: y, x
return x * y
# Область интегрирования:
# x от 0 до 1
# y от 0 до 1-x (треугольная область)
result, error = integrate.dblquad(
f,
0, 1, # пределы по x (нижний, верхний)
lambda x: 0, # нижний предел по y
lambda x: 1 - x # верхний предел по y
)
print("="*50)
print("ПРИМЕР 2: scipy.integrate.dblquad")
print("="*50)
print(f"Функция: f(x, y) = x·y")
print(f"Область: x ∈ [0, 1], y ∈ [0, 1-x]")
print(f"Двойной интеграл ∬ x·y dx dy =
{result:.8f}")
print(f"Оценка ошибки: {error:.2e}")

73. Пример

# Функция: f(x, y, z) = x * y * z
def f(z, y, x): # Порядок аргументов: z, y, x
return x * y * z
# Область интегрирования:
# x от 0 до 1
# y от 0 до 1-x
# z от 0 до 1-x-y
result, error = integrate.tplquad(
f, 0, 1, # пределы по x
lambda x: 0, # нижний предел по y
lambda x: 1 - x, # верхний предел по y
lambda x, y: 0, # нижний предел по z
lambda x, y: 1 - x - y # верхний предел по z
)
print("="*50)
print("ПРИМЕР 3: scipy.integrate.tplquad")
print("="*50)
print(f"Функция: f(x, y, z) = x·y·z")
print(f"Область: x ∈ [0, 1], y ∈ [0, 1-x], z ∈
[0, 1-x-y]")
print(f"Тройной интеграл ∭ x·y·z dx dy dz =
{result:.8f}")
print(f"Оценка ошибки: {error:.2e}")

74. Пример

75. Численное решениe уравнений

76. Численное решениe уравнений

77. Метод деление отрезка пополам (дихотомии)

78. Метод секущих

79. Метод хорд

80. Метод Ньютона (касательных)

81. Метод простых итераций

82. scipy.optimize

83. scipy.optimize

84. scipy.optimize

85. Пример

import numpy as np
import matplotlib.pyplot as plt
import sympy as sp
from scipy import optimize
def bisection(f, a, b, tol=1e-6, max_iter=100):
iter_count = 0
while abs(b - a) >= tol and iter_count < max_iter:
c = (a + b) / 2
if f(c) == 0:
return c, iter_count
elif f(a) * f(c) < 0:
b = c
else:
a = c
iter_count += 1
return (a + b) / 2, iter_count

86. Пример

def secant_method(f, x0, x1, tol=1e-6, max_iter=100):
iter_count = 0
while abs(f(x1)) >= tol and iter_count < max_iter:
x_new = x1 - f(x1) * (x1 - x0) / (f(x1) - f(x0))
x0, x1 = x1, x_new
iter_count += 1
return x_new, iter_count
def chord_method(f, c, x0, tol=1e-6, max_iter=100):
iter_count = 0
while abs(f(x0)) >= tol and iter_count < max_iter:
x0 = x0 - f(x0) * (c - x0) / (f(c) - f(x0))
iter_count += 1
return x0, iter_count

87. Пример

def newton_method(f, df, x0, tol=1e-6, max_iter=100):
iter_count = 0
while abs(f(x0)) >= tol and iter_count < max_iter:
x0 = x0 - f(x0) / df(x0)
iter_count += 1
return x0, iter_count
def simple_iteration(phi, x0, tol=1e-6, max_iter=100):
iter_count = 0
while abs(phi(x0) - x0) >= tol and iter_count < max_iter:
x0 = phi(x0)
iter_count += 1
return x0, iter_count

88. Пример

def draw_function(f, a, b):
x = np.linspace(a, b, 500)
y = f(x)
plt.axvline(x=0, c='gray')
plt.axhline(y=0, c='gray')
plt.plot(x, y, c='steelblue', lw=2)
plt.grid(ls=':')
plt.xlabel('x')
plt.ylabel('f(x)')
plt.title('Поиск уравнения: f(x) = x³ - x - 2')
plt.tight_layout()
plt.show()
def print_result(root_num, root_exact, name_method, iterations='-'):
error = abs(root_num - root_exact)
rel_error = error / abs(root_exact) if root_exact != 0 else error
print(f"| {name_method:<16} | {root_num:^10.8f} | {error:^10.2e} | {rel_error:^10.2%} |
{iterations:^8} |")

89. Пример

x_sym = sp.symbols('x')
expr = x_sym**3 - x_sym - 2
f = sp.lambdify(x_sym, expr, "numpy")
df = sp.lambdify(x_sym, sp.diff(expr, x_sym), "numpy")
# Функция для метода простой итерации: x = phi(x)
# Преобразуем уравнение x^3 - x - 2 = 0 → x = (x+2)^(1/3)
phi = lambda x: (x + 2)**(1/3)
draw_function(f, 0, 2)

90. Пример

print("=" * 70)
print("ЧИСЛЕННОЕ РЕШЕНИЕ НЕЛИНЕЙНЫХ УРАВНЕНИЙ")
print("=" * 70)
print(f"Уравнение: {expr} = 0")
print(f"Интервал локализации [a, b]:")
a = float(input("a = "))
b = float(input("b = "))
x0 = float(input("Начальное приближение: x0 = "))
eps = float(input("Точность: eps = "))
root_exact = float(sp.nsolve(expr, x_sym, x0))
print(f"Точное значение корня: {root_exact:.8f}")

91. Пример

# РУЧНЫЕ РЕАЛИЗАЦИИ МЕТОДОВ
print("=" * 70)
print("РЕЗУЛЬТАТЫ РУЧНОЙ РЕАЛИЗАЦИИ:")
print("=" * 70)
print(f"| {'Метод':<16} | {'Корень':^10} | {'Ошибка':^10} |
{'Отн.ошибка':^10} | {'Итераций':^8} |")
print("-" * 70)
results = {}
try:
root_bisect, iter_bisect = bisection(f, a, b, eps)
results['bisection'] = root_bisect
print_result(root_bisect, root_exact, 'Бисекции', iter_bisect)
except Exception as e:
print(f"| {'Бисекции':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34}
|")

92. Пример

try:
root_secant, iter_secant = secant_method(f, a, b, eps)
results['secant'] = root_secant
print_result(root_secant, root_exact, 'Секущих', iter_secant)
except Exception as e:
print(f"| {'Секущих':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")
try:
root_chord, iter_chord = chord_method(f, a, b, eps)
results['chord'] = root_chord
print_result(root_chord, root_exact, 'Хорд', iter_chord)
except Exception as e:
print(f"| {'Хорд':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")

93. Пример

try:
root_newton, iter_newton = newton_method(f, df, a, eps)
results['newton'] = root_newton
print_result(root_newton, root_exact, 'Ньютона', iter_newton)
except Exception as e:
print(f"| {'Ньютона':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")
try:
root_simple, iter_simple = simple_iteration(phi, x0, eps)
results['simple_iter'] = root_simple
print_result(root_simple, root_exact, 'Простой итерации', iter_simple)
except Exception as e:
print(f"| {'Простой итерации':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34}
|")

94. Пример

# МЕТОДЫ SCIPY.OPTIMIZE
print("=" * 70)
print("ФУНКЦИИ SCIPY.OPTIMIZE:")
print("=" * 70)
# root_scalar с методом brentq (рекомендуемый)
try:
sol_brentq = optimize.root_scalar(f, bracket=[a, b], method='brentq', xtol=eps)
root_brentq = sol_brentq.root
print_result(root_brentq, root_exact, 'brentq (scipy)')
except Exception as e:
print(f"| {'brentq (scipy)':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")
# root_scalar с методом bisect
try:
sol_bisect = optimize.root_scalar(f, bracket=[a, b], method='bisect', xtol=eps)
root_bisect_scipy = sol_bisect.root
print_result(root_bisect_scipy, root_exact, 'bisect (scipy)')
except Exception as e:
print(f"| {'bisect (scipy)':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")

95. Пример

# root_scalar с методом ridder
try:
sol_ridder = optimize.root_scalar(f, bracket=[a, b], method='ridder',
xtol=eps)
root_ridder_scipy = sol_ridder.root
print_result(root_ridder_scipy, root_exact, 'ridder (scipy)')
except Exception as e:
print(f"| {'ridder (scipy)':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")
# root_scalar с методом toms748
try:
sol_toms748 = optimize.root_scalar(f, bracket=[a, b], method='toms748',
xtol=eps)
root_toms748_scipy = sol_toms748.root
print_result(root_toms748_scipy, root_exact, 'toms748 (scipy)')
except Exception as e:
print(f"| {'toms748 (scipy)':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")

96. Пример

# root_scalar с методом newton (требует производную или начальное приближение)
try:
sol_newton_scipy = optimize.root_scalar(f, x0=x0, fprime=df,
method='newton', xtol=eps)
root_newton_scipy = sol_newton_scipy.root
print_result(root_newton_scipy, root_exact, 'newton (scipy)')
except Exception as e:
print(f"| {'newton (scipy)':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")
# newton (упрощённая функция)
try:
root_newton_simple = optimize.newton(f, x0, fprime=df, tol=eps)
print_result(root_newton_simple, root_exact, 'newton (simple)')
except Exception as e:
print(f"| {'newton (simple)':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")

97. Пример

# fsolve (для систем, но можно и для одного уравнения)
try:
root_fsolve = optimize.fsolve(f, x0, xtol=eps)[0]
print_result(root_fsolve, root_exact, 'fsolve (scipy)')
except Exception as e:
print(f"| {'fsolve (scipy)':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")
# root (общий интерфейс)
try:
sol_root = optimize.root(f, x0, method='hybr', tol=eps)
root_root = sol_root.x[0]
print_result(root_root, root_exact, 'root (scipy)')
except Exception as e:
print(f"| {'root (scipy)':<16} | {'ОШИБКА':^10} | {str(e)[:30]:^34} |")
print("=" * 70)

98. Пример

99. Численное решение дифференциальных уравнений

100. Метод Эйлера

Усовершенствованный метод Эйлера
Метод Эйлера-Коши

101. Метод Метод Рунге-Кутты

102. Пример

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import sympy as sp
def euler_method(f, x_span, y0, h):
x = np.arange(x_span[0], x_span[-1] + h, h)
n = len(x)
y = np.zeros(n)
y[0] = y0
for i in range(n - 1):
y[i + 1] = y[i] + h * f(x[i], y[i])
return x, y

103. Пример

def improved_euler_method(f, x_span, y0, h):
x = np.arange(x_span[0], x_span[-1] + h, h)
n = len(x)
y = np.zeros(n)
y[0] = y0
for i in range(n - 1):
y[i + 1] = y[i] + h * f(x[i] + h/2, y[i] + (h/2) * f(x[i], y[i]))
return x, y

104. Пример

def euler_cauchy_method(f, x_span, y0, h):
x = np.arange(x_span[0], x_span[-1] + h, h)
n = len(x)
y = np.zeros(n)
y[0] = y0
for i in range(n - 1):
y_pred = y[i] + h * f(x[i], y[i])
y[i + 1] = y[i] + h/2 * (f(x[i], y[i]) + f(x[i + 1], y_pred))
return x, y

105. Пример

def runge_kutta4_method(f, x_span, y0, h):
x = np.arange(x_span[0], x_span[-1] + h, h)
n = len(x)
y = np.zeros(n)
y[0] = y0
for i in range(n - 1):
k1 = f(x[i], y[i])
k2 = f(x[i] + h/2, y[i] + k1 * h/2)
k3 = f(x[i] + h/2, y[i] + k2 * h/2)
k4 = f(x[i] + h, y[i] + k3 * h)
y[i + 1] = y[i] + h/6 * (k1 + 2*k2 + 2*k3 + k4)
return x, y

106. Пример

def plot_results(results, x_span, h):
x_exact = np.linspace(x_span[0], x_span[-1], 500)
y_exact = analytical_solution(x_exact)
fig, axes = plt.subplots(2, 2, figsize=(14, 8))
axes = axes.flatten()
plt.suptitle(f"Сравнение методов решения ДУ y' = y·x при y(0)=1 и шаге h = {h}",
fontsize=14)
colors = {'Эйлера': 'steelblue', 'Усовершенствованный Эйлера': 'limegreen',
'Эйлера-Коши': 'gold', 'Рунге-Кутта 4': 'violet'}

107. Пример

for idx, (name, (x, y)) in enumerate(results.items()):
ax = axes[idx]
ax.axvline(x=0, c='gray')
ax.axhline(y=0, c='gray')
ax.plot(x, y, 'o-', c=colors[name], lw=2, ms=5, label=f'Метод {name}')
ax.plot(x_exact, y_exact, 'k--', lw=2, label='Точное решение')
y_exact_at_x = analytical_solution(x)
ax.fill_between(x, y, y_exact_at_x, alpha=0.3, color='gray')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_title(f'{name}')
ax.grid(ls=':')
ax.legend(loc=2)
plt.tight_layout()
plt.show()

108. Пример

def f(x, y):
return y * x
# ОСНОВНАЯ ПРОГРАММА
x_span = (0, 2)
y0 = 1.0
h = 0.1
x = sp.symbols('x')
y = sp.Function('y')
de = sp.Eq(y(x).diff(x), y(x)*x)
de_sol = sp.dsolve(de, ics = {y(0): y0})
analytical_solution = sp.lambdify(x, de_sol.rhs, 'numpy')

109. Пример

print("="*80)
print("ЧИСЛЕННОЕ РЕШЕНИЕ ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ")
print("="*80)
print(f"Уравнение: y' = y·x")
print(f"Начальное условие: y(0) = {y0}")
print(f"Аналитическое решение:", de_sol.rhs)
print(f"Интервал: [{x_span[0]}, {x_span[1]}]")
print(f"Шаг интегрирования: h = {h}")
print("="*80)
methods = {
'Эйлера': euler_method,
'Усовершенствованный Эйлера': improved_euler_method,
'Эйлера-Коши': euler_cauchy_method,
'Рунге-Кутта 4': runge_kutta4_method
}

110. Пример

results = {}
for name, method in methods.items():
x, y = method(f, x_span, y0, h)
results[name] = (x, y)
plot_results(results, x_span, h)
y_exact_end = analytical_solution(x_span[-1])
print(f"РЕЗУЛЬТАТЫ ПРИ X = {x_span[-1]}")
print("="*80)
print(f"{'Точное значение':<30} {y_exact_end:^15.8f}")
print("-"*80)
print(f"{'Метод':<30} {'y(x)':^15} {'Ошибка':^15} {'Отн. ошибка':^15}")
print("-"*80)
for name, (x, y) in results.items():
y_end = y[-1]
error = abs(y_end - y_exact_end)
rel_error = error / y_exact_end
print(f"{name:<30} {y_end:^15.8f} {error:^15.2e} {rel_error:^15.2%}")
print("="*80)

111. Пример

112. scipy.integrate: Задача Коши (IVP — Initial Value Problem)

Для решения задачи Коши (ОДУ с начальными условиями) основная функция - solve_ivp.

113. scipy.integrate

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
def ode(t, y):
return y # dy/dt = y
# Решение на [0, 3] с y(0) = 1
sol = solve_ivp(ode, (0, 3), [1],
t_eval=np.linspace(0, 3, 100))
print(sol)
plt.plot(sol.t, sol.y[0])
plt.show()

114. scipy.integrate

y 2 y 4 y x 2 e3 x
y 0 0
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_ivp
import sympy as sp
y0 = [0.0, -5.0]
x_eval = np.linspace(0, 1.5, 1000)
# Аналитическое решение
x = sp.Symbol('x')
y = sp.Function('y')
de = sp.Eq(y(x).diff(x,2)-2*y(x).diff(x)+4*y(x),(x+2)*sp.exp(3*x))
ds = sp.dsolve(de, ics={y(0): y0[0], y(x).diff(x).subs(x, 0): y0[1]})
sp.pprint(ds, wrap_line=False)
y 0 5

115. scipy.integrate

f = sp.lambdify(x, ds.rhs, 'numpy')
plt.figure(figsize=(12, 6))
plt.suptitle(r"Решение ДУ: $y'' - 2y' + 4y = (x+2)e^{3x}$, $y(0)=0,\ y'(0)=-5$",
fontsize=14)
plt.subplot(1, 2, 1)
plt.axvline(x=0, c='gray')
plt.axhline(y=0, c='gray')
plt.plot(x_eval, f(x_eval), c = 'gold', lw=3)
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title('Аналитическое решение: ')
plt.grid(ls=':')

116. scipy.integrate

y 2 y 4 y x 2 e3 x
y y1
y y1 y2
y y2
y2 2 y2 4 y1 x 2 e3 x
y1 y2
3x
y
2
y
4
y
x
2
e
2
1
2
y 0 0
y 0 5
y1 0 0
y2 0 5

117. scipy.integrate

# Численное решение scipy
def ode_system(x, y):
dy0_dx = y[1]
dy1_dx = 2 * y[1] - 4 * y[0] + (x + 2) * np.exp(3 * x)
return [dy0_dx, dy1_dx]
sol = solve_ivp(ode_system, (x_eval[0], x_eval[-1]), y0, t_eval=x_eval, method='RK45')
plt.subplot(1, 2, 2)
plt.axvline(x=0, c='gray')
plt.axhline(y=0, c='gray')
plt.plot(sol.t, sol.y[0], c='steelblue', lw=3)
plt.xlabel('x')
plt.ylabel('y(x)')
plt.grid(ls=':')
plt.title('Численное решение')
plt.tight_layout()
plt.show()

118. scipy.integrate

119. scipy.integrate: Краевая задача (BVP — Boundary Value Problem)

Для решения краевой задачи основная функция - solve_bvp.
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import solve_bvp
def f(x, y):# y'' + y = 0
return np.vstack((y[1], -y[0]))
x_plot = np.linspace(0, np.pi/2, 200)
y_plot = sol.sol(x_plot)[0]
plt.plot(x_plot, y_plot)
# или plt.plot(sol.x, sol.y[0])
plt.show()
def bc(ya, yb): # y(0)=1, y(π/2)=0
return np.array([ya[0] - 1, yb[0]])
x = np.linspace(0, np.pi/2, 30)
y = np.zeros((2, x.size))
sol = solve_bvp(f, bc, x, y)
print(sol)