Добавил:
Upload Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:

Тарасевич Ю.Ю. - Численные методы на Mathcad

.pdf
Скачиваний:
108
Добавлен:
16.03.2015
Размер:
796.56 Кб
Скачать

f2( 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

TOL := 10− 3

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