Численные методы

В этом листке мы познакомимся с разными методами численного дифференцирования и интегрирования, а также сравним их. Но перед этим нам потребуется научиться разбирать математические функции, записанные текстом, строить синтаксические деревья, вычислять значения такого выражения и символьно вычислять производную. :)

Вспоминаем, как рисовать графики

A: Вспомнить matplotlib

Для разминки нарисуйте график функции $$2\cos(2x) + 1/2 \cdot \sin(x^2)$$ на отрезке от -10 до 10 по 401 точке. Вспоминать, что это, можно здесь.

Если бы точек было 41, а не 401, то график бы выглядел так:
img

B: Вычислить выражение при помощи eval

В питоне есть отличная встроенная функция eval, которая позволяет исполнить переданный код и вернуть значение получившегося выражения. Например, вычислить $2\sin(\pi/6)$ можно и так:

from math import sin, pi
x = pi / 6
print(eval('sin(x) * 2'))

В отличие о явного программирования, строку 'sin(x)*2' можно получить снаружи и пока не парсить.

Итак, на вход в первой сточке дана запись математической функции, в которой кроме арифметики, чисел, скобок, переменной $x$ может участвовать также функции $\sin$, $\cos$, $\exp$, $\log$ и другие математические функции из модуля math. Во второй строчке дано действительное число $x$. Вычислите значение функции в этой точке.

sin(1/2*exp(x**0.5)) + cos(x) 1.79
0.7270445075667562

C: Нарисовать график при помощи eval

На вход в первой сточке дана запись математической функции от одной переменной $x$.
Во второй строчке даны два действительных числа — границы отрезка, на котором нужно нарисовать график функции.
В третьей сточке дано целое число — количество точек.
Нарисуйте график данной функции на данном отрезке.

sin(1/x) -17.9 17.9 40
img

Обратная польская запись

В математике существует древняя традиция помещать оператор между операндами $(x+y)$, а не после операндов $(xy+)$. Форма с оператором между операндами называется инфиксной записью. Форма с оператором после операндов называется постфиксной, или обратной польской записью в честь польского логика Я. Лукасевича (1958), который изучал свойства этой записи.

Примеры инфиксной и постфиксной записи:

ИнсфикснаяПостфиксная
$3 + 4 \times 2 / (1 - 5)^2$3 4 2 × 1 5 - 2 ^ / +
$\exp(-1/2\times x)$1 ~ 2 / x × exp, здесь ~ — это унарный минус
$ 3 + 4 \times 2 / (1 - 5 ) ^ {2 ^ 3}$3 4 2 × 1 5 - 2 3 ^ ^ / +
$\sin ( \max ( 2, 3 ) / 3 \times \pi )$2 3 max 3 / 𝜋 × sin

Обратная польская запись имеет ряд преимуществ перед инфиксной записью при выражении алгебраических формул. Во-первых, любая формула может быть выражена без скобок. Во-вторых, она удобна для вычисления. В-третьих, инфиксные операторы имеют приоритеты, которые произвольны и нежелательны. Например, мы знаем, что $ab+c$ значит $(ab)+c$, а не $a(b+c)$, поскольку произвольно было определено, что умножение имеет приоритет над сложением. Но имеет ли приоритет сдвиг влево перед логической операцией И? Кто знает? Обратная польская запись позволяет устранить такие недоразумения.

Алгоритм сортирующей станции

Для построения обратной польской записи будем использовать алгоритм сортирующей станции, идея которого предложена Э. Дейкстра (E.W. Dijkstra). Посмотрите его оригинальный текст 1961 года, это забавно.

x

На рисунке схематично показана железная дорога из Нью-Йорка в Калифорнию с ответвлением, ведущим в Техас. Каждый операнд/оператор формулы представлен одним вагоном. Поезд движется на запад (налево). Перед стрелкой каждый вагон должен останавливаться и узнавать, должен ли он двигаться прямо в Калифорнию, или ему нужно будет по пути заехать в Техас. Вагоны, содержащие числа и переменные, всегда направляются в Калифорнию и никогда не едут в Техас. Вагоны, содержащие все прочие операторы, должны перед прохождением стрелки «оценивать ситуацию». Иногда в результате этой «оценки» сначала потребуется перегнать несколько вагонов из Техаса в Калифорнию, а ещё парочку вагонов сбросить с рельсов в пропасть :)

В определении того, что должен делать вагон, основная тонкость. Но с ней достаточно легко разобраться, добавляя возможные операторы небольшими порциями.

Вычисление выражения в обратной польской записи

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

Вот пример вычисления значения выражения $(1 + 2)\times 4 + 3$, который в ОПЗ имеет вид 1 2 + 4 × 3 +.
x

D1: ОПЗ — плюс-минус

Дано арифметическое выражение с целыми числами и арифметическими операциями +, - (без скобок и унарных плюсов/минусов). Все операторы/операнды разделены пробелами. Преобразуйте его в обратную польскую запись.

3 + 4
3 4 +
1 - 2 + 3 + 4 - 5
1 2 - 3 + 4 + 5 -

D2: ВЫЧ — плюс-минус

Дана запись выражения в обратной польской записи с целыми числами и арифметическими операциями +, - (без унарных плюсов/минусов). Вычислите его значение

3 4 +
7
1 2 - 3 + 4 + 5 -
1

E1: ОПЗ — умножение-деление

Дано арифметическое выражение с целыми числами и арифметическими операциями +, -, *, / (без скобок и унарных плюсов/минусов). Все операторы/операнды разделены пробелами. Преобразуйте его в обратную польскую запись.

2 + 2 * 2
2 2 2 * +
1 + 2 * 3 * 4 + 6 / 2 + 1
1 2 3 * 4 * + 6 2 / + 1 +

E2: ВЫЧ — умножение-деление

Дана запись выражения в обратной польской записи с целыми числами и арифметическими операциями +, -, *, / (без унарных плюсов/минусов). Вычислите его значение

2 2 2 * +
6
1 2 3 * 4 * + 6 2 / + 1 +
29.0

F1: ОПЗ — скобки

Дано арифметическое выражение с целыми числами и арифметическими операциями +, -, *, /, в котором могут быть скобки (но без унарных плюсов/минусов). Все операторы/операнды разделены пробелами. Преобразуйте его в обратную польскую запись.

(2 + 2) * 2
2 2 + 2 *
(1 + 2) * (3 + 4) / (5 + 6 / 2)
1 2 + 3 4 + * 5 6 2 / + /

G1: ОПЗ — унарные плюс-минус

Дано арифметическое выражение с целыми числами и арифметическими операциями +, -, *, /, в котором могут быть скобки и унарные плюсы/минусы. Обратите внимание, что последовательно унарных плюсов/минусов вычисляется права налево. Все операторы/операнды разделены пробелами. Преобразуйте его в обратную польскую запись. Для обозначения унарного минуса пишите u-, для унарного плюса — u+.

- 2 + (- 2 + 2) - - - + 2
2 u- 2 u- 2 + + 2 u+ u- u- -

G2: ВЫЧ — унарные плюс-минус

Дана запись выражения в обратной польской записи с целыми числами и арифметическими операциями +, -, *, /, u-, u+ Вычислите его значение

2 u- 2 u- 2 + + 2 u+ u- u- -
-4

H1: ОПЗ — возведение в степень

Дано выражение с числами и операциями +, -, *, /, ^. Обратите внимание, несколько возведений в степень вычисляются справа налево. Преобразуйте его в обратную польскую запись.

2 + 2 * 2 ^ 2
2 2 2 2 ^ * +
1 ^ 2 ^ 3 ^ 4
1 2 3 4 ^ ^ ^

H2: ВЫЧ — возведение в степень

Дана запись выражения в обратной польской записи с числами и арифметическими операциями +, -, *, /, u-, u+, ^ Вычислите его значение

2 2 2 2 ^ * +
10

I1: ОПЗ — функции с одним параметром

Дано выражение с числами, операциями +, -, *, /, ^, sin, cos, exp, ... и прочими математическими функциями из библиотеки math. Гарантируется, что имя функции — это более одной латинской буквы. Корректность функции проверять не нужно. Преобразуйте его в обратную польскую запись.

sin(asin(2 - 1)) * atan(tan(1))
2 1 - asin sin 1 tan atan *

I2: ВЫЧ — функции с одним параметром

Дана запись выражения в обратной польской записи с числами и арифметическими операциями +, -, *, /, u-, u+, ^ и функциями exp ln sqrt acos asin atan cos sin tan с одним параметром. Вычислите его значение

2 1 - asin sin 1 tan atan *
1.0

J1: ОПЗ — переменные

Дано выражение с числами, операциями +, -, *, /, ^, sin, cos, exp, ... и, возможно, переменными. Гарантируется, что имя функции — это более одной латинской буквы, а переменная — просто одна буква (не обязательно латинская, например, 𝜋). Преобразуйте его в обратную польскую запись.

x * y + sin(2 * 𝜋)
x y * 2 𝜋 * sin +

J2: ВЫЧ — переменные

Дана запись выражения в обратной польской записи с переменными. В следующей строчке дано количество переменных $v$. Затем идёт $v$ пар строчек — имя переменной и её значение. Вычислите значение выражения с подставленными константами.

x y * 2 𝜋 * sin + 3 x 11 y 22 𝜋 3.141592
241.9999986928204

K1★: ОПЗ — без пробелов

Дано выражение с числами, операциями +, -, *, /, ^, sin, cos, exp, ... и, возможно, переменными. Гарантируется, что имя функции — это более одной латинской буквы, а переменная — просто одна буква (не обязательно латинская, например, 𝜋). Однако между операндами может не быть пробелов вовсе или может быть больше одного пробела. Преобразуйте его в обратную польскую запись.

x*y+sin( 2 * 𝜋 )
x y * 2 𝜋 * sin +

L1★: ОПЗ — пропущены умножения

Дано выражение с числами, операциями +, -, *, /, ^, sin, cos, exp, ... и, возможно, переменными. В некоторых случаях оператор умножения * может быть опущен. Например, $2x$, $(2+x)(2-x)$, $2\sin(x)$, $2\pi r$. Гарантируется, что имя функции — это одна из следующих строк: exp ln log sqrt acos asin atan cos sin tan. Всё остальное — это переменные, причём переменные никогда не стоят «впритык» к имени функции.

x y sin(2𝜋)cos((1+x)(2+2x y))
x y * 2 𝜋 * sin 1 x + 2 2 x * y * + * cos *

Синтаксические деревья для математических формул

После того, как мы научились строить обратную польскую запись и вычислять по ней значения, уже достаточно легко строить синтаксические деревья. Например, для формулы $$ x y \sin(2 \pi) \cos((1+x)(2+2xy)) $$ обратная польская запись имеет вид x y * 2 𝜋 * sin 1 x + 2 2 x * y * + * cos *, синтаксическое дерево выглядит так:

Абстрактное синтаксическое дерево — в информатике ориентированное дерево, в котором внутренние вершины сопоставлены с операторами, а листья — с операндами (то есть с переменными и константами). В синтаксическом дереве для математических формул в листьях записаны константы и переменные, в вершинах дерева записаны функции, и в каждую внутреннюю вершину входит ровно столько рёбер, сколько параметров имеет данная функция. Скажем в бинарные операции входит две стрелки, в тригонометрические функции и унарные операции — по одной.

Синтаксические деревья также позволяют эффективно вычислять значения выражения при помощи достаточно простой рекурсивной процедуры. Однако деревья позволяют легко выполнять более сложные операции. Например, можно упрощать деревья, выполняя операции, в которых участвуют только константы. Или использовать дистрибутивность. Нам же в этом листке синтаксические деревья помогут символьно вычислять производную функции.

Построение синтаксического дерева

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

class Node:
    def __init__(self, op, parms):
        self.op = op
        self.parms = parms

nx = Node(None, ['x'])
ny = Node(None, ['y'])
nxy = Node('*', [nx, ny])
n2 = Node(None, ['2'])
n𝜋 = Node(None, ['𝜋'])
n2𝜋 = Node('*', [n2, n𝜋])
nsin2𝜋 = Node('sin', [n2𝜋])
...

Будем хранить в атрибуте op оператор или None, если это константа. А в списке parms будем хранить либо список из единственной константы, либо список из нод, которые являются параметрами данного оператора.

Строить синтаксическое дерево можно, вычисляя значение обратной польской записи. Только в стек нужно класть не результаты вычислений, а экземпляры класса Node.

M1: Строим синтаксическое дерево по ОПЗ

Дано выражение с обратной польской записью математической формулы. Постройте по нему синтаксическое дерево. Соберите все его рёбра в список строк вида ребёнок родитель. Отсортируйте его как список строк и выведите.

1 x + 2 y + * sin
* sin + * + * 1 + 2 + x + y +

M2: Вычисляем значение в синтаксическом дереве

Дано синтаксическое дерево и значения переменных. Нужно вычислить значение выражения.

В первой строчке дано количество $N$ вершин дерева. Затем идут $N$ строк с описаниями вершины в таком формате:
номер_вершины|оператор|константа_или_список_номеров_вершин_параметров_через_запятую
В следующей строчке дано количество переменных $v$. Затем идёт $v$ пар строчек — имя переменной и её значение. Вычислите значение выражения с подставленными константами.

8 0||1 1||x 2|+|0,1 3||2 4||y 5|+|3,4 6|*|2,5 7|sin|6 2 x 2.141592653589793 y 0.5
1.0
Это выражение равно $sin(\pi \cdot 2.5) = \sin(\pi/2) = 1$.

M3: Собираем формулу из дерева

Дано синтаксическое дерево. Нужно вывести его в человеко-читаемом виде.

В первой строчке дано количество $N$ вершин дерева. Затем идут $N$ строк с описаниями вершины в таком формате:
номер_вершины|оператор|константа_или_список_номеров_вершин_параметров_через_запятую
Сборку нужно выполнить при помощи рекурсивной процедуры. У функций обязательно писать скобки.

8 0||1 1||x 2|+|0,1 3||2 4||y 5|+|3,4 6|*|2,5 7|sin|6
sin((1+x)*(2+y))

Таблица производных

Вычисление производной — важнейшая операция в дифференциальном исчислении. Ниже список формул для нахождения производных от некоторых функций.

Константы и степени

$$(c)' = 0$$ $$(x)' = 1$$ $$(x^c)' = cx^{c-1}, c\ne0$$

Экспоненты и логарифмы

$$(e^x)'=e^x$$ $$(c^x)'=c^x \cdot \ln c$$ $$(\ln x)' = \frac{1}{x}$$ $$(\log_a x)' = \frac{1}{x \ln a}$$

Тригонометрия

$$(\sin x)' = \cos x$$ $$(\cos x)' = -\sin x$$

Произведение, отношение и композиция

$$(fg)' = f' g + fg'$$ $$\left(\frac{f}{g}\right)' = \frac{f' g - f g'}{g^2}$$ $$(f(g(x))' = f'(g(x)) \cdot g'(x)$$

Символьное вычисление производной

Если у нас уже построено синтаксическое дерево для математической функции, то по нему достаточно легко посчитать производную. Это можно сделать рекурсивно, используя формулы выше.

N: Производная многочлена

Дано выражение, содержащее арифметику и переменную $x$. Вычислите символьно производную этой функции. Оборачивайте все не-константы в скобки. Раскрывать скобки, убирать скобки, убирать умножения на 1 не нужно. Прибавление нулей и умножения на нули нужно убирать.

1 + 2 * x
(2)*(1)
(1 + x) * (1 + 2*x)
(1)*(1+2*x)+(1+x)*(2)
(1 + x * x) / (3 * x - 1)
((2*x)*(3*x-1)-(1+x*x)*(3*1))/(3*x-1)^2

O: Производная произвольного выражения

Дано выражение, содержащее арифметику, функции exp, ln, sqrt, acos, asin, atan, cos, sin, tan, и переменную $x$. Вычислите символьно производную этой функции. Оборачивайте все не-константы в скобки. Раскрывать скобки, убирать скобки, убирать умножения на 1 не нужно. Прибавление нулей и умножения на нули нужно убирать.

Ответ будет проверяться не по совпадению текста, а по совпадению значения функции в наборе точек.

sin(x)
cos(x)
(1 + x) * (1 + 2*x)
(1)*(1+2*x)+(1+x)*(2)
sin(2 * x ** 2) * sqrt(log(x))
2*2*cos(2*x**2)*sqrt(log(x))+sin(2*x**2)*1/(2*x*sqrt(log(x)))

Символьное вычисление производной и вычисление значение в библиотеке sympy

Ниже пример всех шагов:

from sympy.parsing.sympy_parser import parse_expr
from sympy import diff
from sympy.abc import x

expr = 'sin(2 * x ** 2) * sqrt(log(x))'
f = parse_expr(expr)
df = diff(f, x)

x_0 = 179
print(f)
print(df)
print(f.subs({x: x_0}).evalf())
print(df.subs({x: x_0}).evalf())

Соглашение для графиков

Для упрощения изложения мы всегда будет рисовать графики только на отрезке $[-10; 10]$ по $N$ точкам. Число $N$ всегда будет первым параметром теста.

P: График функции и её точной производной

На вход в первой строке даётся число $N$ — количество точек. Во второй строке даётся функция от переменной $x$. Нарисуйте график этой функции и её производной.

51 sin(2 * sqrt(x + 11))
ing

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

Теперь представим себе, что нам неизвестна сама функция, а только известны её значения $f_i$ в $N$ точках вида $x_k = a + k\cdot h$, где $k=0,\ldots,N-1$. Тогда может потребоваться как-нибудь оценить её значения в промежуточных точках. Эта процедура называется интерполя́ция. Точки $x_i$ называются узлами, а их совокупность — сеткой. Число $h$ называется шагом сетки.

Про интерполяцию есть много разной красивой математики, но в рамках этого листка мы будем использовать только линейную интерполяцию. То есть будем заменять настоящую неизвестную функцию на $(N-1)$-звенную ломаную, проходящую через точки $(x_i, f_i)$. Тогда значение в точке $x\in[x_i, x_{i+1}]$ получается равным $$ f(x) = f_i + \frac{f_{i+1}-f_i}{x_{i+1}-x_{i}} (x - x_i). $$

Когда вы рисуете график функции в matplotlib, вы как раз автоматом получаете вместо точного графика линейную интерполяцию.

Q: Интерполяция

На вход в первой строке даётся число $N$ — количество узлов сетки (на отрезке $[-10, 10]$). Во второй строке даётся функция от переменной $x$. В третьей точке дано количество дополнительных узлов $K$ для линейной интерполяции. То есть на каждом из $(N-1)$ отрезков нужно равномерно добавить по $K$ узлов, значения функции в которых нужно вычислить при помощи линейной интерполяции.

Необходимо вывести $N + K\cdot(N-1)$ чисел — значения функции в узлах сетки и линейная интерполяция в промежуточных.

3 (x+5)*(x+5) 3
25 25.0 25.0 25.0 25 75.0 125.0 175.0 225

Численное дифференцирование

Нам неизвестна сама функция, а только известны её значения $f_i$ в $N$ точках вида $x_k = a + k\cdot h$, где $k=0,\ldots,N-1$. Самый наивный способ численно посчитать производную — метод односторонней разности: в узлах сетки положим $$ df_i = \frac{f_{i+1} - f_{i}}{x_{i+1}-x_{i}} $$ А в промежуточных точках вычислим при помощи линейной интерполяции. Значение в самой правой точке придётся взять равным значению в предпоследней точке.

R: Метод односторонней разности

На вход в первой строке даётся число $N$ — количество узлов сетки (на отрезке $[-10, 10]$). Во второй строке даётся функция от переменной $x$. В третьей точке дано количество дополнительных узлов $K$ для линейной интерполяции. Вычислите численно производную при помощи метода односторонней разности. В промежуточных узлах используйте линейную интерполяцию.

Постройте два графика на отрезке $[-10, 10]$ по $N + K\cdot(N-1)$ точкам: численное значение производной и точное.

3 (x+5)*(x+5) 3
ing
15 sin(2 * sqrt(x + 11)) 63
ing

Численное дифференцирование: метод двусторонней разности

Оказывается, что гораздо лучшее значение производной даёт — метод двусторонней разности: в узлах сетки положим $$ df_i = \frac{f_{i+1} - f_{i-1}}{x_{i+1}-x_{i-1}} $$ Значение в концевых точках придётся посчитать при помощи метода односторонней разности (правосторонней и левосторонней). А в промежуточных точках вычислим при помощи линейной интерполяции.

S: Метод двусторонней разности

На вход в первой строке даётся число $N$ — количество узлов сетки (на отрезке $[-10, 10]$). Во второй строке даётся функция от переменной $x$. В третьей точке дано количество дополнительных узлов $K$ для линейной интерполяции. Вычислите численно производную при помощи метода односторонней разности. В промежуточных узлах используйте линейную интерполяцию.

Постройте три графика на отрезке $[-10, 10]$ по $N + K\cdot(N-1)$ точкам:
численное значение производной при помощи метода односторонней разности;
численное значение производной при помощи метода двусторонней разности;
значение символьно посчитанной производной.

Чтобы изменить толщину линии добавьте методу ax.plot параметр linewidth. Чтобы создать легенду добавьте методу ax.plot параметр label, а также добавьте вызов метода ax.legend().

15 sin(4 * sqrt(x + 11)) 63
ing

T: Метод двусторонней разности — 2

Всё то же самое, что в предыдущей задаче. Но вместо самой производной постройте два графика с разностью численно вычисленной и символьно вычисленной.

Чтобы изменить толщину линии добавьте методу ax.plot параметр linewidth. Чтобы создать легенду добавьте методу ax.plot параметр label, а также добавьте вызов метода ax.legend().

41 sin(4 * sqrt(x + 11)) 63
ing

Численное дифференцирование: магия 5 точек

Оказывается, что можно считать ещё лучше по четырём точкам, если расстояние $h$ между соседними узлами одинаковое. В узлах сетки положим $$ df_i = \frac{f_{i-2} - 8 f_{i-1} + 8 f_{i+1} - f_{i+2}}{12(x_{i+1}-x_{i}} $$ Откуда взялись эти магические константы 1, -8, 8, -1? Тут потребуется матан и ряды Тейлора. Пока поверьте на слово :)

U: Производная по 4-м точкам

Всё то же самое, что в предыдущей задаче, только добавляется третий график. Ещё отрежьте два крайних приближения, там будут самые большие ошибки. И да, выводите модули :)

Постройте три графика на отрезке $[-10 + \frac{40}{N-1}, 10 - \frac{40}{N-1}]$ по $(N-4) + K\cdot(N-5)$ точкам: численное значение производной и точное.

41 sin(4 * sqrt(x + 11)) 63
ing