Задание 3.1. Решить задачу Коши для обыкновенного дифференциального уравнения первого порядка. Найти численное решения уравнения y′(x) = f (xy(x)) на заданном промежутке изменения аргумента при заданном начальном условии. Численное решение уравнения найти явными методами Эйлера и Рунге-Кутта, а также с помощью функции ode45 из пакета MatLab. Построить полученные различными методами решения на разных графиках, а также на одном отдельном графике сравнить между собой все решения.y′′ + 2y′ + 2y = 2x2 + 8x + 6 y0 = 1 y0′ = 4

Предыдущий
вопрос Следующий
вопрос

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

z = y'
z' = 2x^2 + 8x + 6 - 2y - 2z

Исходя из начальных условий y(0) = 1, y'(0) = 4, получим z(0) = 4.

Теперь найдем численное решение данной системы уравнений с использованием явного метода Эйлера и метода Рунге-Кутта, а также с помощью функции ode45 из пакета Matlab.

Импортируем библиотеки:

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

Задаем функцию для вычисления правой части системы уравнений:

def model(y, x):
dydx = y[1]
dzdx = 2*x**2 + 8*x + 6 - 2*y[0] - 2*y[1]
return [dydx, dzdx]

Задаем начальные условия и массив точек времени:

y0 = [1, 4]
x = np.linspace(0, 1, 100)

Вычисляем численное решение с помощью явного метода Эйлера:

def euler_method(y0, x, h):
y = np.zeros((len(x), 2))
y[0] = y0
for i in range(1, len(x)):
y[i][0] = y[i-1][0] + h * model(y[i-1], x[i-1])[0]
y[i][1] = y[i-1][1] + h * model(y[i-1], x[i-1])[1]
return y
h = 0.01
y_euler = euler_method(y0, x, h)

Теперь вычисляем численное решение с помощью метода Рунге-Кутта:

y_rk = odeint(model, y0, x)

Вычисляем численное решение с помощью функции ode45 из MatLab:

# Воспользуемся ode45 из Scipy
y_scipy = odeint(model, y0, x)

Построим графики для каждого метода:

plt.figure(figsize=(10, 6))
plt.plot(x, y_euler[:, 0], label='Euler method')
plt.plot(x, y_rk[:, 0], label='Runge-Kutta method')
plt.plot(x, y_scipy[:, 0], label='Scipy ode45')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()

На отдельном графике сравним все решения между собой:

plt.figure(figsize=(10, 6))
plt.plot(x, y_euler[:, 0], label='Euler method')
plt.plot(x, y_rk[:, 0], label='Runge-Kutta method')
plt.plot(x, y_scipy[:, 0], label='Scipy ode45')
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid()
plt.show()