Тарасевич Ю.Ю. - Численные методы на Mathcad
.pdff2( z) yi
z,xi
f3( z) |
yi |
z,xi |
3.3.3. Аппроксимация линейной комбинацией функций |
Mathcad предоставляет пользователям встроенную функцию linfit для аппроксимации данных по методу наименьших квадратов линейной комби- нацией произвольных функций.
31
Функция linfit имеет три аргумента:
∙ вектор x – x–координаты заданных точек, ∙ вектор y – y–координаты заданных точек,
∙ функция F – содержит набор функций, который будет использоваться
для построения линейной комбинации. |
|
|||||||
Задаем |
функцию F |
(аппроксимирующая функция ищется |
в виде |
|||||
1 |
|
|
|
|
1 |
|
|
|
b·x |
2 |
: |
F( x ) |
x |
1 |
|
|
|
a x +1 + |
|
x2 |
|
|
|
|||
|
|
|
|
|
|
|
|
|
S linfit ( x , y , F ) |
|
|
|
|
||||
0.802 |
|
|
|
|
|
|
||
S = |
|
|
|
|
|
|
|
|
0.176 |
|
|
|
|
|
|
||
Определяем аппроксимирующую функцию: f4( x ) F( x ).S |
|
|||||||
Вычисляем дисперсию: |
|
|
|
|||||
S4 |
( n |
1 . |
f4 xi |
yi |
2 |
|
||
|
|
2) |
|
|
|
|
||
|
|
|
|
i |
|
|
S4 = |
0.975 |
f4( z) |
|
|
|
|
|
|
|
|
yi |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
z,xi |
|
32 |
|
|
|
|
|
|
|
|
3.3.4. Аппроксимация функцией произвольного вида
Теперь построим аппроксимирующую функцию дробно–рационального
типа f(x) = |
ax2 |
. Для этого воспользуемся функцией genfit. Функция |
|
b + x |
|||
|
|
имеет следующие параметры:
∙x, y – векторы, содержащие координаты заданных точек,
∙F – функция, задающая искомую функциональную n–параметрическую зависимость и частные производные этой зависимости по параметрам.
∙v – вектор, задающий начальные приближения для поиска параметров.
|
|
|
|
|
|
|
|
|
|
|
|
|
u |
.z2 |
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
0 |
|
|
|
|
|
|
|
|
|
|
|||||
|
|
|
|
|
|
|
|
|
|
|
|
|
u1 |
|
|
|
|
|
z |
|
|
|
||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||
F( z , u) |
|
|
|
|
|
|
z2 |
|
|
|
|
|||||||||||||
|
|
|
|
|
|
u1 |
|
|
|
|
|
z |
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
u |
|
.z2 |
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
0 |
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
u |
|
|
|
|
|
z 2 |
|
|
|
||
|
|
|
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
|
|
|
|||||
1 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
v |
|
|
|
|
|
15 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
2.146 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
S |
|
genfit( x , y , v, F ) |
S = |
|||||||||||||||||||||
|
||||||||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
20.85 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
Поскольку нулевой элемент функции F содержит искомую функцию, определяем функцию следующим образом: f5( z ) F( z , S )0
Вычисляем среднее квадратичное отклонение
S5 |
1 . |
f5 xi |
yi |
2 |
S5 = 0.581 |
( n |
2) |
i |
|
|
|
|
|
|
|
|
33
f5( z) |
|
|
|
|
|
|
|
|
|
|
|
|
|
yi |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
z,xi |
|
|
|
|
|
|
|
|
Функция genfit имеется не во всех реализациях Mathcad'а. Возможно, |
|||||||||||||
однако, решить задачу, проведя линеаризацию. |
|
|
|
|
|
|
|
|
|||||
Заданная функциональная зависимость может быть линеаризована вве- |
|||||||||||||
дением переменных |
z = 1 |
и t = 1 . Тогда z = |
1 |
+ b |
t . |
|
|
|
|||||
|
|
|
|
y |
x |
|
|
a |
|
a |
|
|
|
Определим матрицы коэффициентов нормальной системы (см. книгу |
|||||||||||||
[8] из списка литературы) |
|
|
|
|
|
|
|
|
|
||||
x |
i |
3.y |
i |
|
x |
i |
4 |
|
x |
i |
2.y |
i |
|
|
|
|
|
|
|
|
|
|
|||||
i |
|
|
|
|
i |
|
|
|
|
i |
|
|
|
e |
|
|
2 |
o |
|
2 |
|
|
|
|
2 |
|
|
x . |
|
y |
x |
i |
.y |
i |
|
y |
|
||||
i |
|
|
i |
|
|
|
|
|
|
i |
|
i |
i |
i |
Находим коэффициенты функции, решая систему матричным методом,
|
|
o 1.e |
|
|
1.218 |
|
|
|
|||
d |
|
d = |
|||
|
|||||
|
|||||
|
|
|
|
|
15.517 |
|
|
|
|
Определяем функцию: f5( z )
z2.d0
d1 z
Вычислим стандартное отклонение
34
S5 |
1 . |
f5 xi |
yi |
2 |
S5 = 0.827 |
( n |
2) |
i |
|
|
|
|
|
|
|
|
Обратите внимание! Мы получили другие коэффициенты! Вспомни- те, задача на нахождение минимума нелинейной функции, особенно не- скольких переменных, может иметь несколько решений.
Стандартное отклонение больше, чем в случае аппроксимации полино- мами, поэтому следует остановить свой выбор на аппроксимации полино- мом.
Представим результаты аппроксимации на графиках
f5( z) |
yi |
z,xi |
В тех случаях, когда функциональная зависимость оказывается доста- |
точно сложной, может оказаться, что самый простой способ нахождения |
коэффициентов – минимизация функционала Ф "в лоб". |
35
Глава 4. Вычисление определенных интегралов
4.1. Метод Ромберга
Пусть требуется вычислить определенный интеграл на интервале [a;b].
b
ò f (x)dx
a
Далеко не всегда задача может быть решена аналитически. В частности, численное решение требуется в том случае, когда подынтегральная функ- ция задана таблично. Для численного интегрирования подынтегральную функцию аппроксимируют какой-либо более простой функцией, интеграл от которой может быть вычислен. Обычно в качестве аппроксимирующей функции используют полином. В случае полинома нулевой степени метод численного интегрирования называют методом прямоугольников, в случае полинома первой степени – методом трапеций, в случае полинома второй степени – методом Симпсона. Все эти методы являются частными случая-
ми квадратурных формул Ньютона-Котеса.
Итак, в методе трапеций подынтегральную функцию аппроксимируют полиномом первой степени, то есть прямой линией. Это значит, что вместо площади криволинейной трапеции мы будем искать площадь прямоуголь- ной трапеции. Приближенное значение интеграла равно
b |
f (b) - f (a) |
|
||||||
ò f (x)dx » |
(b - a) |
|||||||
|
||||||||
a |
2 |
|
|
|
|
|||
Погрешность этой формулы равна O(h |
3 |
f |
¢¢ |
|||||
|
) . |
|||||||
Обозначим R(1;1) = |
h |
( f (a) - f (b)) , где |
|
h = b - a . Смысл введенного |
||||
|
|
|||||||
2 |
|
|
|
|
|
|
||
обозначения станет ясен несколько позже. |
|
|
|
|
Оценку значения интеграла можно сделать более точной, если разбить интервал на n частей и применить формулу трапеций для каждого такого
интервала
æ 1 |
ö |
||
hç |
|
( f (a) + f (b)) + f (a + h) + f (a + 2h) +L+ f (a + (n -1)h))÷ . |
|
2 |
|||
è |
ø |
Если разбить интервал на две части, то есть уменьшит шаг в два раза h = (b - a)2 , то оценка для величины интеграла будет иметь вид
|
h |
( f (a) + f (b)) + hf (a + h) = |
1 |
2n −1 |
|
R(2;1) = |
R(1;1) + hå f (a + ih) |
||||
|
2 |
||||
2 |
|
i=1 |
36
В данном случае суммирование включает только один элемент. Обра- тите внимание, в новую оценку вошла старая оценка. Нам потребовалось определять значение функции только в новых узлах.
Если имеется 2n подынтервалов, то
|
|
|
|
|
|
|
h |
|
|
|
|
|
|
2n −1 |
||||
|
|
|
R(n +1;1) = |
( f (a) + f (b))+ hå f (a + ih) |
||||||||||||||
2 |
||||||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
|
b − a |
|
i=1 |
||||
|
|
|
|
|
|
|
|
|
|
h = |
|
|
|
|||||
Если n=0, то |
|
|
2n |
|||||||||||||||
|
|
|
|
|
|
|
||||||||||||
|
|
|
|
|
h |
|
|
|
|
|
||||||||
|
|
|
|
|
|
R(1;1) = |
|
( f (a) + f (b)) |
||||||||||
2 |
||||||||||||||||||
Если n=1, то |
|
|
|
|
|
|
|
|
|
|||||||||
|
|
|
|
|
|
|
|
|
|
|
||||||||
|
R(2;1) = |
h |
( f (a) + f (b)) |
+ hf (a + h) = |
1 |
R(1;1) + hf (a + h) |
||||||||||||
|
|
|||||||||||||||||
2 |
|
|
|
|
|
2 |
|
|||||||||||
Если n=2, то |
|
|
|
|
|
|
|
|
|
|
|
|||||||
R(3;1) = |
h |
( f (a) + f (b))+ hf (a + h) + hf (a + 2h) + hf (a + 3h) = |
||||||||||||||||
|
||||||||||||||||||
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||||
= |
1 |
R(2;1) + hf (a + h) + hf (a + 3h) |
||||||||||||||||
|
||||||||||||||||||
2 |
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
|
||
Вообще, справедливо рекуррентное соотношение |
||||||||||||||||||
|
|
|
|
|
|
|
|
1 |
|
|
|
|
2n −1 |
|||||
|
|
|
|
R(n +1;1) = |
|
R(n;1) + hå f (a + (2i −1)h) |
||||||||||||
|
|
|
|
2 |
||||||||||||||
|
|
|
|
|
|
|
|
|
|
|
i=1 |
|||||||
Полученное соотношение называют рекурсивной формулой трапеций и |
часто применяют для вычисления определенных интегралов. Преимущест- во этой формулы состоит в том, что при увеличении числа подынтервалов функцию нужно вычислять только во вновь добавленных точках. К сожале- нию, с помощью этой формулы нельзя получить сколь угодно точное зна- чение интеграла. Во-первых, при увеличении числа разбиений объем вы- числений стремительно возрастает; во-вторых, на каждом шаге накаплива- ется ошибка округлений. Для дальнейшего уточнения значения интеграла можно сделать следующий шаг – экстраполировать полученную последо-
вательность значений на случай бесконечного числа точек или что то же самое, на случай нулевого шага. Такой подход называется методом Ром-
берга.
Метод Ромберга заключается в том, что полученные оценки значения интеграла экстраполируют на случай бесконечного числа разбиений (вели-
чины шага равной нулю) по рекуррентной формуле |
|
R(n +1, m +1) = R(n +1, m) + R(n +1,m) − R(n, m) |
(4.1) |
4m −1
То есть строится следующий треугольник
R(1,1)
37
R(2,1) R(2,2) R(3,1) R(3,2) R(3,3)
R(4,1) R(4,2) R(4,3) R(4,4)
R(5,1) R(5,2) R(5,3) R(5,4) R(5,5) ,
в котором первый столбец состоит из значений интеграла, полученных при последовательном удвоении числа интервалов. Второй столбец – ре-
зультат уточнения значений первого столбца по рекуррентной формуле (4.1). Третий столбец – уточненные значения интеграла на основе второго столбца и т.д.
Формула (4.1) может быть получена различными способами. Можно, например, воспользоваться методом Невиля. Пусть имеется набор точек (xi , yi ) . Обозначим Pi полином нулевой степени, проходящий через i-ю
точку. Обозначим Pi(i+1) полином первой степени, проходящий через точки i и i+1. Совершенно аналогично будет означать P12K(n−1)n полином n–1 степе- ни, проходящий через все n точек. Легко убедиться, что
Pi(i +1)K(i+m) = |
(x − xi+m )Pi(i+1)K(i+m−1) + (xi |
− x)P(i+1)(i+2)K(i+m) |
||||
|
b − a |
xi − xi+m |
|
|
||
|
|
|
|
|||
В нашем случае x = |
. В качестве y |
выступают R(i;m) . Мы хотим |
||||
2i−1 |
||||||
|
i |
i |
|
|
||
|
|
|
|
|
получить значение интеграла в пределе h → 0 , поэтому x = 0 .
4.2 Вычисление определенных интегралов
Для вычисления определенного интеграла необходимо выбрать знак интеграла из палитры или набрать его нажатием клавиши &. После этого следует вписать пределы интегрирования, подынтегральную функцию и переменную интегрирования. Mathcad успешно справляется с большинст- вом интегралов, в том числе, с несобственными. Точность вычислений ре- гулируется встроенной переменной TOL. По умолчанию ее значение уста-
новлено . Ниже приводится несколько примеров успешного вы- числения несобственных интегралов, интеграла от быстро осциллирующей функции и интеграла от ступенчатой функции.
ó |
− 1 |
|
|
|
|
|
|
ô |
|
1 |
dx = 1 |
ó |
∞ |
|
2 |
|
|
− x |
|||||
ô |
2 |
|
ô |
x× e |
dx = 0.5 |
||
õ |
|
x |
|
õ |
|
||
− ∞ |
|
0 |
|
|
|||
ó |
20 |
|
(x3) dx = 1.822´ 10− 3 |
ó |
1 |
|
|
ô |
sin |
ô |
F(x) dx = 1 |
||||
õ10 |
|
|
õ− 1 |
|
|
38
ì0, если x<0
Здесь F(x) = í
î1, если x ³ 0
Зависимость результата от заданной точности вычислений представле-
на ниже
ó∞ |
|
1 |
|
|
ô |
|
|
dx = 3.14159265369356 |
|
|
|
|
||
ô |
1 |
2 |
|
|
õ− ∞ |
+ x |
|
||
|
|
|
|
|
TOL:= 10− 6 |
|
|||
ó∞ |
|
1 |
|
|
ô |
|
|
dx = 3.14159265358979 |
|
|
|
|
||
ô |
1 |
2 |
|
|
õ− ∞ |
+ x |
|
||
|
|
|
|
Для этого примера результат может быть получен также в символьном виде. Для этого вместо знака равенства необходимо нажать знак символи- ческого равенства Ctrl+.
ó∞ |
|
1 |
|
ô |
|
dx ® p |
|
|
|
||
ô |
1 |
2 |
|
õ− ∞ |
+ x |
|
|
|
|
|
В то же время в некоторых случаях несобственные интегралы вычис- ляются неправильно.
ó |
1 |
1 |
dx = 1.376´ 103 |
ô |
|
||
|
|
||
ô |
2 |
|
|
|
|
x |
|
õ− 1 |
|||
Хотя очевидно, что |
|||
ó |
1 |
1 |
|
ô |
|
dx ® ¥ |
|
|
|
||
ô |
2 |
|
õ− 1 x
39
Глава 5. Решение дифференциальных уравнений
5.1. Обыкновенные дифференциальные уравнения
Введение
Пусть необходимо найти решение уравнения |
|
|
y′ = |
f (x, y) |
(5.1) |
с начальным условием y(x0 ) = |
y0 . Такая задача называется задачей |
|
Коши. Разложим искомую функцию |
y(x) в ряд вблизи точки |
x0 и ограни- |
чимся первыми двумя членами разложения y(x) = y(x0 ) + y′(x)(x − x0 ) +K .
Учтя |
уравнение |
(5.1) |
и |
обозначив |
x − x0 = h , |
получаем |
y(x) = y(x0 ) + f (x0 , y0 ) |
x. Эту формулу можно |
применять |
многократно, |
|||
находя значения функции во все новых и новых точках. |
|
|||||
|
|
|
yi+1 = yi |
+ f (xi , yi )h |
|
(5.2) |
Такой метод решения обыкновенных дифференциальных уравнений назы- вается методом Эйлера. Геометрически метод Эйлера означает, что на ка- ждом шаге мы аппроксимируем решение (интегральную кривую) отрезком касательной, проведенной к графику решения в начале интервала. Точность метода невелика и имеет порядок h. Говорят, что метод Эйлера – метод первого порядка, то есть его точность растет линейно с уменьшением шага h.
Существуют различные модификации метода Эйлера, позволяющие увеличить его точность. Все они основаны на том, что производную, вы- численную в начале интервала, заменяют на среднее значение производной на данном интервале. Среднее значение производной можно получить (ко- нечно же только приближенно) различными способами. Можно, например, оценить значение производной в середине интервала [xi , xi +1 ] и использо-
вать его для аппроксимации решения на всем интервале
x |
|
|
= x + |
|
h |
|
|
|
|
|
|
1 |
|
|
|
|
|
|
|
|
|||
i+ |
|
i |
2 |
|
|
|
|
|
|
||
2 |
|
|
|
|
|
|
|
|
|||
y |
|
|
= y + |
h |
f (x , y |
) |
|
|
|||
1 |
|
|
|
|
|||||||
i+ |
|
i |
2 |
i |
i |
|
|
|
|||
2 |
|
|
|
|
|
|
|
||||
y |
|
|
= y + hf (x |
, y |
1 |
) |
|||||
i+1 |
i |
|
|
|
i + 1 |
|
i + |
|
|||
|
|
|
|
|
|
|
2 |
|
|
2 |
|
Можно также оценить среднее значение производной на интервале y%i+1 = yi + hf (xi , yi )
40