Инженерные расчеты на Python: библиотеки и примеры

Решение системы нелинейных уравнений

Как найти пересечение двух кривых, заданных неявно?

import numpy as np
from scipy.optimize import fsolve

# Система: x^2 + y^2 = 4, x*y = 1
def system(vars):
    x, y = vars
    eq1 = x**2 + y**2 - 4
    eq2 = x*y - 1
    return [eq1, eq2]

initial = [1, 1]
solution = fsolve(system, initial)
print(f"Решение: x = {solution[0]:.6f}, y = {solution[1]:.6f}")
print(f"Проверка: x^2+y^2 = {solution[0]**2+solution[1]**2:.6f}, x*y = {solution[0]*solution[1]:.6f}")
  

Решение: x = 0.517638, y = 1.931852
Проверка: x^2+y^2 = 4.000000, x*y = 1.000000
  

Пояснение: fsolve ожидает функцию, возвращающую список невязок. Начальное приближение подбирается исходя из физического смысла задачи. При неудачном выборе начальной точки алгоритм может найти другое решение или не сойтись. Рекомендуется визуализировать функции для выбора хорошего старта.

Подбор параметров модели по экспериментальным данным (curve_fit)

Как аппроксимировать набор точек экспоненциальной функцией с оценкой погрешности?

import numpy as np
from scipy.optimize import curve_fit
import matplotlib.pyplot as plt

# Модельные данные: y = A * exp(-B*x) + шум
xdata = np.linspace(0, 4, 50)
y_true = 2.5 * np.exp(-1.3 * xdata)
y_noise = 0.2 * np.random.normal(size=xdata.size)
ydata = y_true + y_noise

# Определяем модель
def model(x, A, B):
    return A * np.exp(-B * x)

# Выполняем подгонку
popt, pcov = curve_fit(model, xdata, ydata, p0=[1, 1])
print(f"A = {popt[0]:.3f}, B = {popt[1]:.3f}")
print("Ковариационная матрица:")
print(pcov)

# Визуализация
plt.plot(xdata, ydata, 'b.', label='Экспериментальные точки')
plt.plot(xdata, model(xdata, *popt), 'r-', label='Аппроксимация')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.show()
  

A = 2.487, B = 1.289
Ковариационная матрица:
[[0.0045 0.0012]
 [0.0012 0.0007]]
(график не выводится в консоли, но в среде Jupyter отобразится)
  

Численное интегрирование с контролем точности

Как вычислить определенный интеграл с заданной точностью, если аналитический первообразной нет?

from scipy.integrate import quad
import numpy as np

# Интеграл от sin(x)/x на [0, 10]
def f(x):
    return np.sin(x) / x

result, error = quad(f, 0, 10, epsabs=1e-8, epsrel=1e-6)
print(f"Результат: {result:.10f}, оценка погрешности: {error:.3e}")
  

Результат: 1.6583475943, оценка погрешности: 1.084e-13
  

Пояснение: quad использует адаптивный алгоритм Симпсона или Кронрода. Параметры epsabs и epsrel задают требуемую абсолютную и относительную погрешность. Для несобственных интегралов можно указать границы с бесконечностью (например, 0, np.inf). Проблема может возникнуть, если подынтегральная функция имеет особенности на концах интервала тогда нужно задать точки разрыва через параметр points.

Обработка сигналов: быстрое преобразование Фурье

Как выделить основные частоты в зашумленном сигнале?

import numpy as np
import matplotlib.pyplot as plt
from scipy.fft import fft, fftfreq

# Параметры сигнала
fs = 1000  # частота дискретизации, Гц
T = 1.0    # длительность, с
N = int(fs * T)
t = np.linspace(0, T, N, endpoint=False)

# Сумма двух синусоид + шум
f1, f2 = 50, 120
signal = 0.6 * np.sin(2*np.pi*f1*t) + 0.4 * np.sin(2*np.pi*f2*t)
noise = 0.2 * np.random.normal(size=N)
total = signal + noise

# БПФ
yf = fft(total)
xf = fftfreq(N, 1/fs)[:N//2]
amp = 2.0/N * np.abs(yf[:N//2])

# Поиск пиков (простейший пороговый детектор)
threshold = 0.1
peaks = xf[amp > threshold]
print("Частоты с амплитудой выше порога:", peaks)

# Визуализация
plt.plot(xf, amp)
plt.xlabel('Частота, Гц')
plt.ylabel('Амплитуда')
plt.grid()
plt.show()
  

Частоты с амплитудой выше порога: [ 50. 120.]
(график спектра)