Практикум j практическое примщенше численных методов

1 т
= /(*»+>-»^п+1)-
В частности, при т = 3 имеем схему
l
/(*»+!, У^П+1),

lyⁿ⁺¹ - 18у^п + 9у^п~¹ - 2у^п~² _
6т
которая имеет третий порядок аппроксимации.

3. 
   Упражнения
   

Упражнение 12.1 Напишите программу для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений явным ме­тодом Рунге—Кутта четвертого порядка. Продемонстрируйте работоспо­собность этой программы при решении задачи Коши
s
d²u dt²
in(w), 0 < t < 47Г,
, ч du , _ч
u(0) = l, — (0) = 0.
В модуле rungeKutta функция rungeKuttaO реализует решение задачи Коши для системы ОДУ методом Рунге—Кутта четвертого порядка.
Модуль rungeKutta^ import numpy as np
def rungeKuttaCf, tO, yO, tEnd, tau):
и и и
Solve the initial value problem y’ = f(t,y) by 4th-order Runge-Kutta method. t0,y0 are the initial conditions, tEnd is the terminal value of t, tau is the step.
и и и
def increment(f, t, y, tau): kO = tau * f(t,y)
kl = tau * f(t + tau/2., у + kO/2.)
k2 = tau * f(t + tau/2., у + kl/2.)
k3 = tau * f(t + tau, у + k2)
return (kO + 2.*kl + 2.*k2 + k3) / 6. t = [] у = []
t.append(tO) у.append(yO) while tO < tEnd:
tau = min(tau, tEnd - tO)
yO = yO + increment(f, tO, yO, tau)
tO = tO + tau t.append(tO) у.append(yO)
return np.array(t), np.array(y)
П

dy₂
dt


sin^x), 0 < t < 47г.

ри приближенном решении модельной задачи Коши для уравнения второго порядка сначала переходим от одного уравнения второго порядка к системе из двух уравнений
Р



Рис. 12.1 Приближенное решение задачи при т = 0.25

ешение задачи при заданном шаге интегрирования т обеспечивает следую­щая программа.
import numpy as np import math as mt import matplotlib.pyplot as pit from rungeKutta import rungeKutta def f(t, y):
f = np.zeros((2),’float ’)
f СО] = у[1] f[l] = - mt. sin(y [0]) return f tO = 0.
tEnd = 4.*np.pi
yO = np.array([l., 0.])
tau =0.25
t, у = rungeKutta(f, tO, yO, tEnd, tau) for n in range(0, 2): r = y[:, n] st = ’$y$’ sg =
if n == 1:
st = ^,$\\frac{d y}{dt}$’ sg =
pit.plot(t, r, sg, label=st) pit.Iegend(loc=0) pit.xlabel(’$t$’) pit.grid(True) pit.show()
В рассматриваемой задаче (колебания маятника) при малых у имеем sin (у) « у и решение уравнение имеет период равный 2тг. Нелинейность проявляется, в частности, в увеличении периода (см. рис. 12.1).
Упражнение 12.2 Напишите программу для численного решения задачи Коши для системы обыкновенных дифференциальных уравнений с использо­ванием двухслойной схемы с весом при решении системы нелинейных урав­нений на новом временном слое методом Ньютона. Используйте эту про­грамму для решения задачи Коши (модель Лотка—Вольтерра)
^ = 2/1 - 2/12/2, = “2/2 + 2/12/2, 0 < t < 10,
2/i (0) = 2, г/2 (0) = 2.
Реализация двухслойной (одношаговой) схемы с весом проводится функцией oneStepO в модуле oneStep. Для решения системы нелинейных уравнений используется модуль newton.
Модуль оneStep:н
import numpy as np
from newton import newton
def oneStep(f, tO, yO, tEnd, nTime, theta):
и и и
Solve the initial value problem y* = f(t,y) by one-step methods implisit method.
t0,y0 are the initial conditions, tEnd is the terminal value of t, nTime is the number of steps.
и и и
tau = (tEnd - tO) / nTime def fl(yl):
fl = yl - tau * theta * f(t0+tau,yl) \
- yO - tau * (1.-theta) * f(t0+tau,y0) return fl t = []
у = []
t.append(tO) у.append(yO) for i in range(nTime):
r = yO - tau * f(t0+tau,y0) yl, iter = newton(fl, r) yO = yl tO = tO + tau t.append(tO) у.append(yO)
return np.array(t), np.array(y)
Решение модельной задачи Коши при постоянном шаге по времени обеспечи­вает следующая программа.
import numpy as np import matplotlib.pyplot as pit from oneStep import oneStep def f(t, y):
f = np.zeros((2),’float’) f [0] = у [0] - у [0] *y [1] f [1] = -y [1] + y[0]*y[l] return f tO = 0. tEnd = 10.
yO = np.array([2., 2.]) nTime = 50 theta =0.5
t, у = oneStep(f, tO, yO, tEnd, nTime, theta) for n in range(0, 2): r = у [:, n]
st = ’$y_l$’
sg =
if n == 1:


Рис. 12.2 Решение задачи при в = 0.5