Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Кафедра прикладной математики

Творческая работа

по предмету «Численные методы»

на тему

Численные методы решения уравнений

2011

Содержание

Решение алгебраических, нелинейных и трансцендентных уравнений

Метод половинного деления (дихотомия)

Метод простых итераций

Метод касательных (Ньютона)

Метод секущих

Численные методы вычисления определённых интегралов

Метод левых прямоугольников

Метод правых прямоугольников

Метод средних прямоугольников

Метод трапеций

Приближение функций

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

Аппроксимация

Степенная регрессия

Показательная регрессия

Численные методы решения дифференциальных уравнений

Постановка задачи

Метод Эйлера

Общая формулировка методов Рунге-Кутты

Обсуждение методов порядка 4

«Оптимальные» формулы

Условия порядков для методов Рунге-Кутты

Оценка погрешности и сходимость методов Рунге-Кутты

Строгие оценки погрешности

Главный член погрешности

Оценка глобальной погрешности

Оптимальный выбор шага

Список использованной литературы

Решение алгебраических, нелинейных и трансцендентных уравнений

Пусть задана непрерывная функция f(х) и требуется найти все или некоторые корни уравнения

f(x)=0. (1)

Эта задача распадается на несколько задач. Во-первых, надо исследовать количество, характер и расположение корней. Во-вторых, найти приближенные значения корней. В-третьих, выбрать из них интересующие нас корни и вычислить их с требуемой точностью.

Первая и вторая задачи решаются аналитическими и графическими методами.

Когда ищутся только действительные корни уравнения, то полезно составить таблицу значений f(x). Если в двух соседних узлах таблицы функция имеет разные знаки, то между этими узлами лежит нечетное число корней уравнения (по меньшей мере, один). Если эти узлы близки, то, скорее всего, корень между ними только один. Но выявить по таблице корни чётной кратности*^* Кратными корнями называются совпадающие по значению корни. Например, известно, что уравнение y=x² имеет корень x=0, но правильнее утверждать, что это уравнение имеет 2 кратных корня x₁=0 и x₂=0. сложно. По таблице можно построить график функции у=f(х) и графически найти точки его пересечения с осью абсцисс. Этот способ более нагляден и дает неплохие приближенные значения корней. Во многих задачах техники такая точность уже достаточна. В технике еще популярны графические методы решения уравнений (номография). Построение графика позволяет выявить даже корни чётной кратности.

Иногда удается заменить уравнение (1) эквивалентным ему уравнением (х)=(х), в котором функции y1=(х) и y2=(х) имеют несложные графики. Например, уравнение хsinх--1=0 удобно преобразовать к виду sinx=l/x. Абсциссы точек пересечения этих графиков будут корнями исходного уравнения.

Приближенные значения корней уточняют различными итерационными методами. Рассмотрим наиболее эффективные из них.

Метод половинного деления (дихотомия)

Пусть мы нашли такие точки a и b что f(a)f(b)0, т. е. на отрезке [a,b] лежит не менее одного корня уравнения. Найдем середину отрезка xc=(a+b)/2 и вычислим f(xc). Из двух половин отрезка выберем ту, для которой f(xc)f(a или b)0, т.е. отрезок на котором функция меняет знак. Затем новый отрезок опять делим пополам и выберем ту половину, на концах которой функция имеет разные знаки, и т. д. (рис. 1).

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Если требуется найти корень с точностью , то продолжаем деление пополам до тех пор, пока длина отрезка не станет меньше 2. Тогда середина последнего отрезка даст значение корня с требуемой точностью. Дихотомия проста и очень надежна: к простому корню она сходится для любых непрерывных функций f(x), в том числе недифференцируемых; при этом она устойчива к ошибкам округления. Скорость сходимости невелика: за одну итерацию точность увеличивается примерно вдвое, т. е. уточнение трех цифр требует 10 итераций (т.к. длина отрезка, на котором лежит корень, после 10 итераций равна 1/210=1/102410-3). Зато точность ответа гарантируется.

Перечислим недостатки метода.

1. Для начала расчета надо найти отрезок, на котором функция меняет знак.

2. Если в этом отрезке несколько корней, то заранее неизвестно, к какому из них сойдется процесс (хотя к одному из них сойдется).

3. Метод неприменим к корням четной кратности.

4. Для корней высокой нечетной кратности он сходится, но менее точен и хуже устойчив к ошибкам округления, возникающим при вычислении f(x).

5. Наконец, на системы уравнений дихотомия не обобщается.

Утверждение 1. С помощью данного метода невозможно найти корни чётной кратности.

Доказательство.

Чётно кратный корень это корень уравнения вида

(x+a)2n=0, где n - целое, n[0,]. (2)

Решением этого уравнения будет корень x=-a кратности 2n. В общем виде уравнение может иметь как чётно, так и нечётно кратные корни. Можно записать общий вид уравнения имеющего (k+m) только действительных корней так:

(x+x1)2n1(x+x2)2n2…(x+xk)2nk(x+xk+1)2n(k+1)+1(x+xk+2)2n(k+2)+1…(x+xk+m)2n(k+m)+1=0, (3)

где n1,…,n(k+m) [0,] - целые числа; x1 x2… xk+m.

В уравнении (3) k чётно кратных и m нечётно кратных корней. Оно раскладывается на (k+m) уравнений, из которых легко получаются корни. Если задать начальный отрезок [-x1-r,-x1+r], где r - мало, и проверить условие смены знака функции на его границах, то обнаружим, что знак не меняется в силу чётности степени. А если аналогично проверить нечётно кратные корни, то получим обратную ситуацию.

Следствие 1.

Если корень имеет чётную кратность, то на границах бесконечно малого отрезка с центром в этом корне функция имеет одинаковые знаки.

Следствие 2.

Если корень имеет нечётную кратность, то на границах бесконечно малого отрезка с центром в этом корне функция имеет разные знаки.

Пусть на заданном отрезке [a,b] лежит 1 корень чётной кратности, тогда в силу следствия 1 на границах отрезка знак меняться не будет, что означает остановку выполнения итераций и недостижение необходимой точности. Если же на отрезке [a,b] лежит 1 чётно кратный корень и 1 нечётно кратный корень, то чётно кратный корень будет просто игнорирован методом, т.к. условие смены знака являющееся также основным условием, с помощью которого определяется корень на текущем полуотрезке, в силу следствия 1 не выполнится. Следовательно, чётно кратный корень не может быть найден с помощью данного метода.

Утверждение 2. Если на концах начального отрезка значения функции имеют один знак, то метод может не сойтись, то есть, возможно, ни один из корней не будет найден с заданной точностью.

Доказательство.

Первым вариантом постоянства знака функции на границах отрезка является отсутствие корня на нём, поэтому исключим этот случай как тривиальный, будем считать, что на отрезке хотя бы один корень существует. Вторым вариантом - существование чётного количества корней.

Если f(a)f(b)0, то продолжать итерации невозможно, т.к. условие смены знака не подтверждается. Если же, тем не менее, на первом шаге не проверять условие смены знака и разделить отрезок пополам, то может возникнуть ситуация, в которой корни распределяться по чётному количеству в каждой половине отрезка. А чётное количество корней означает чётное количество пересечений оси Ox, даже если существуют кратные корни.

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

Утверждение 3. Если на концах начального отрезка значения функции имеют разные знаки, то будет найден с заданной точностью один из корней лежащих на нём.

Доказательство.

В силу утверждения 1 будем рассматривать только корни нечётной кратности. Так как функция меняет знак на концах отрезка, предположим, f(a)?0, f(b)?0. Тогда если f(xc)?0, то для дальнейшего приближения выберем отрезок [xc,b], т.к. f(b)f(xc)?0. Если же f(xc)?0, то для дальнейшего приближения выберем отрезок [a,xc], т.к. f(a)f(xc)?0.Для второго случая, когда f(a)?0, f(b)?0 аналогично доказывается существование одного из полуотрезков, на котором функция меняет знак. Из чего следует, что после каждой итерации для одного из полуотрезков условие смены знака обязательно будет выполнено. Следовательно, нет причин для остановки итерационного процесса, который завершится лишь по достижении заданной точности.

Дихотомия применяется тогда, когда требуется высокая надежность счета, а скорость сходимости малосущественна.

Метод простых итераций

Основной принцип метода заключается в том, что уравнение (1) представляется в виде:

x=ц(x), (4)

где ц(x) можно определить многими способами, например, так:

ц(x)=x-бf(x), б=const, или

ц(x)=x+ш(x)f(x),

где ш(x) - произвольная, непрерывная, знакопостоянная функция на отрезке [a,b].

Метод простых итераций в силу (4) определяется следующей рекурсивной формулой:

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

xn+1=ц(xn), где n=0,1,2,… (5)

Здесь n имеет смысл номера итерации, x0 - некоторое начальное приближение. Из (5) видно, что если xn>xт, то этот предел и есть корень уравнения (рис. 2).

Пусть в окрестности точки xт (xт-Д,xт+Д), где Д>0 функция ц(x) удовлетворяет условию Липшица:

|ц(x2)-ц(x1)|?q|x2-x1| (6)

для любых x2,x1(xт-Д,xт+Д),

0<q<1, (7)

при этом x0(xт-Д,xт+Д), (8)

причём, (xт-Д,xт+Д) [a,b].

В связи с этим допущением можно сделать несколько утверждений.

Утверждение 1.

Полученные с помощью (5) xn(xт-Д,xт+Д) для любого целого n?0.

Доказательство

В силу (8)

|x0-xт|<Д,

из (5),(6) и (7) получим

|x1-xт|=|ц(x0)-ц(xт)|?q|x0-xт|,

|x2-xт|=|ц(x1)-ц(xт)|?q|x1-xт|?q2| x0-xт|,

|xn-xт|=|ц(xn-1)-ц(xт)|?q|xn-1-xт|?q2| xn-2-xт|?…? qn| x0-xт|<Д, (9)

Из последнего неравенства |xn-xт|<Д следует, что для любого целого n?0 верно утверждение 1.

Утверждение 2.

Последовательность {xn} сходится при n>? к пределу xт, являющемуся корнем уравнения.

Доказательство.

В силу (5) и (6):

|xn+m-xn|=|ц(xn+m-1)-ц(xn-1)|?q|xn+m-1-xn-1|?q2|xn+m-2-xn-2|?…?qn|xm-x0| (10)

Следовательно, в силу (7) для любого целого

m?0

Значит в силу признака Коши*^* Признак Коши. Для того чтобы последовательность {x_n} имела предел, необходимо и достаточно, чтобы для любого как угодно малого положительного числа е существовал такой номер N=N_е, при котором для любых n?N и всех целых положительных m выполнялось неравенство |x_n+m-x_n|<е, другими словами предел существует . Докажем теперь, что xlim=xт. Из (10) следует (для определённости возьмём m=1), что , т.е. при n>? xlim=xn+1=xn, а в силу (5) xn+1=ц(xn) => xlim=xn=ц(xn) => xlim=ц(xlim). Из этого равенства следует, что xlim есть корень уравнения (4), т.е. xlim=xт.

Утверждение 3.

На интервале (xт-Д,xт+Д) существует только 1 корень уравнения (2).

Доказательство.

Пусть существует 2 корня

x1 и x2?x1, тогда в силу (4) x1=ц(x1), x2=ц(x2).

Из (6) следует |x2-x1|=|ц(x2)-ц(x1)|?q|x2-x1|, т.е. |x2-x1|?q|x2-x1|.

Из этого следует, что q=1, а это противоречит (7). Следовательно, корень только 1.

Следует так же отметить, что если ц(x) имеет на интервале (xт-Д,xт+Д) непрерывную производную, то (6) выполняется, когда

|ц?(x)|?q => |ц?(x)|<1 для x(xт-Д,xт+Д). (11)

Если |ц?(xт)|<1, x0(xт-Д,xт+Д), то итерации сойдутся.

Если |ц?(xт)|>1, то в силу непрерывности производной |ц?(x)|>1 и на некотором интервале (xт-Д1,xт+Д1), поэтому итерации не сойдутся к корню. Если же |ц?(xт)|<1, но x0(xт-Д,xт+Д), т.е. |ц?(x0)|>1, то в данном случае о сходимости процесса можно судить только после дополнительного анализа функции ц(x).

Как и при поиске решения методом дихотомии будем считать задачу выполненной, если найденное некоторое значение xч[xт-е,xт+е], где е - заданная точность. Для определения того, когда можно прекратить итерации, т.е. когда достигнута заданная точность, подробнее рассмотрим неравенство (9). По сути, нам необходимо добиться выполнения следующего неравенства:

|xn+1-xт| ? е, (12)

а в силу того, что из (9) можно получить неравенство

|xn+1-xт| ? qn+1|x0-xт|,

где q можно определить как для любого целого n?1, выведем условие достижения заданной точности (12). Введём обозначения

дn+1=|xn+1-xт|. Из (9) ,

а так же очевидно, что д0-дn+1=|x0-xn+1|=о, тогда:

д0=дn+1+о

дn+1=gд=g(дn+1+о)



. (13)

Неравенство (13) является условием остановки процесса итераций, т.е. условием достижения заданной точности.

Метод касательных (Ньютона)

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Метод Ньютона называют также методом касательных и методом линеаризации. Суть метода заключается в том, что в точке приближения к функции строится касательная (Рис. 3). Следующая точка приближения - это точка пересечения полученной прямой с осью Ox. Процесс продолжается вплоть до достижения заданной точности.

Из рисунка очень легко получить итерационную формулу метода, используя геометрический смысл производной. Если f(x) имеет непрерывную производную f'(x)?0, тогда получим



Аналогично получаем x2, x3, и т.д. Таким образом, можем записать общую формулу:

Метод Ньютона можно рассматривать, как частный случай метода простых итераций, если задать

.

В этом случае

Условие сходимости метода простых итераций можно переписать для метода Ньютона следующим образом:

(14)

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

Метод секущих

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

В данном методе, в отличие от метода Ньютона, проводятся не касательные, а секущие (Рис. 4). Из рисунка легко получить итерационную формулу:

. (15)

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

В знаменателе формулы (15) стоит разность значений функции. Вдали от корня это не существенно; но вблизи корня, особенно корня высокой кратности, значения функции малы и очень близки. Возникает потеря значащих цифр, приводящая к «разболтке» счёта. Это ограничивает точность, с которой можно найти корень; для простых корней это ограничение не велико, а для кратных может быть существенным.

От «разболтки» страхуются так называемым приёмом Гаврика. Выбирают не очень малое е, ведут итерации до выполнения условия |xn+1-xn|<е и затем продолжают расчёты до тех пор, пока |xn+1-xn| убывают. Первое же возрастание обычно означает начало «разболтки»; тогда расчёт прекращают и последнюю итерацию не используют. Все ограничения по сходимости итераций для данного метода такие же, как и в методах простых итераций и Ньютона. А вот определение достижения заданной точности, как видно из описания метода, затруднительно, и, даже, возможна ситуация, когда из-за «разболтки» счёта заданная точность не будет достигнута никогда. При использовании метода секущих в явном виде определить точность трудно, поэтому используют косвенный метод. Считают, что вблизи корня |xn+1-xn||xт-xn+1|. Конечно эта оценка весьма примерна, но при больших n (в идеале при n>?) это так и есть. Построим блок-схему алгоритма данного метода с использованием приёма Гаврика.

Численные методы вычисления определённых интегралов

Метод левых прямоугольников

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Эти методы численного интегрирования основываются на геометрическом смысле интеграла. Как известно интеграл есть площадь криволинейной трапеции ограниченной подынтегральной функцией f(x) на отрезке [a,b]. Для вычисления интеграла I отрезок [a,b] разбивают на n отрезков длиной h. На каждом отрезке криволинейную трапецию приближают прямоугольником, так как его площадь можно легко вычислить. Затем суммируют все полученные площади, получая тем самым приближённое значение интеграла. На рис. 5 проиллюстрирован, так называемый, метод левых прямоугольников. Его название объясняется тем, что высота прямоугольника f(x) вычисляется в левой границе отрезка h. Выведем формулу для приближённого вычисления интеграла I.

Как видно из рисунка площадь первого слева прямоугольника S0=f(x0)h. Площадь следующего S1=f(x1)h. Легко заметить, что площадь i-го прямоугольника Si=f(xi)h. Всего таких прямоугольников n, нумерация их ведётся от 0 до n-1. Таким образом, приближённое значение интеграла, полученное этим методом, вычисляется по формуле:

. (16)

Оценим погрешность формулы (16) на отрезке [xi,xi+h]. Для этого разложим функцию f(x) по формуле Тейлора, выбирая xi за центр разложения и предполагая наличие у функции требуемых по ходу рассуждений непрерывных производных.

f(x)=f(xi)+(x-xi)f'(xi)+… (17)

В формуле (17) для определённости отбросим члены, содержащие производные высших порядков, т.е. 2-ю и выше. Тогда получим:

, т.е.

(18)

Чтобы получить общую погрешность метода на отрезке [a,b] просуммируем все RЛi:

(19)

Поскольку в (17) были отброшены члены, содержащие более высокие степени длины интервала, то выражение остаточного члена (19) является асимптотическим, т.е. выполняющимся при h>0 с точностью до членов более высокого порядка малости. Но для справедливости этой оценки необходимо существование непрерывной f'(x); если f'(x) кусочно-непрерывна, то удаётся сделать лишь мажорантную оценку

(20)

Эту формулу можно записать иначе

, (21)

где n - количество отрезков, на которые разбит отрезок [a,b].

Перед началом вычислений по данному методу необходимо определиться с h или, что, по сути, то же самое, с n. Пусть известны границы отрезка [a,b], задана подынтегральная функция f(x) и точность е, с которой должен быть вычислен интеграл. Чтобы определить h и n воспользуемся простым и очевидным неравенством:

|RЛ|?е, (22)

откуда получаем для непрерывной f'(x):

(23)

либо если учесть, что получим

. (24)

Для кусочно-непрерывной f'(x) с учётом (20) и (21) справедливо следующее

(25)

(26)

Метод правых прямоугольников

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Этот метод похож на предыдущиу. Отличие в том, что высота прямоугольников вычисляется по правой границе (Рис. 6).

Выводы формул для данного метода аналогичны предыдущему. Основные отличия заключаются в нумерации. Формула метода правых прямоугольников выглядит следующим образом:

. (27)

Рассуждения и конечные формулы, связанные с определением погрешности метода и выбором ширины основания прямоугольников аналогичны тому, что получено в предыдущем методе.

Метод средних прямоугольников

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Чтобы уменьшить погрешность методов левых и правых прямоугольников был предложен метод средних, т.е. метод в котором высота прямоугольника вычисляется в середине отрезка h (Рис. 7). Обращаясь к рисунку легко увидеть, что площади прямоугольников вычисляются по следующим формулам:

(28)

Оценим погрешность формулы (28) на отрезке [xi,xi+h]. Прежде всего, получим разложение по формуле Тейлора для данного метода. Центр разложения в данном случае будет точка . Тогда получим:

(29)

В формуле (29) отбросим члены, содержащие 3-ю и выше производные. Получаем:

Ii-IСр.i=RСр.i

, т.е.

. (30)

Аналогично тому, как были получены формулы (20), (21), (23), (24), (25) можно получить формулы, по которым оцениваются n и h, для данного метода. Запишем окончательные формулы для непрерывной f''(x).

(31)

(32)

Для кусочно-непрерывной f''(x) получаем мажорантные оценки:

(33)

(34)

Сравнивая формулы (19) и (30) погрешностей методов приходим к выводу, что погрешность метода средних во много раз ниже погрешности метода левых или правых прямоугольников, т.е. метод средних во много раз точнее и для достижения заданной точности требует меньше машинного времени. Поэтому блок-схему построим именно для метода средних.

Введём обозначение. [n] - целая часть n, полученная путём отбрасывания дробно части. Для вычисления интеграла дано е, [a,b], f(x), f''(x). Прежде всего, через е необходимо получить n. Пользоваться формулой (32) мы не будем, т.к. в ней необходимо вычислять интеграл. Воспользуемся мажорантной оценкой (34). Для этого перед вычислением интеграла необходимо найти , но это легко осуществимо с помощью небольшого циклического процесса. Для него понадобиться знание количества разбиений , на которых будем вычислять значения второй производной. Обозначим это количество n1.

Метод трапеций

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Метод трапеций основан на том, что криволинейная трапеция приближается прямолинейной (Рис. 8). Т.е. площади вычисляются по следующей формуле:

Таким образом, получаем общую формулу трапеций:

. (35)

Теперь оценим погрешность метода. Вывод формулы погрешности аналогичен выводу (30), поэтому приведём сразу окончательную оценку погрешности для непрерывной f''(x):

; (36)

для кусочно-непрерывной f''(x):

(37)

Оценивая n и h, получим:

(38)

(39)

Для кусочно-непрерывной f''(x):

(40)

 (41)

Как видно из полученных формул погрешность метода средних примерно вдвое меньше погрешности метода трапеций. Очевидно, что если функция определена на всём интервале, лучше пользоваться методом средних, метод трапеций используют обычно для функций определённых только в узлах сетки. Знаки главного члена погрешности у формул трапеций и средних разные. Поэтому, если есть расчёты по обеим формулам, то точное значение интервала лежит, как правило, в вилке между ними. Деление этой вилки как 2:1 даёт уточнённый результат близкий к тому, который получается при использовании более точного метода Симпсона.

Приближение функций

Как известно, функцию можно задать многими способами. Самые распространённые из них - это аналитический, графический и табличный. В первых двух способах функция обычно определена на бесконечном числе точек, третий же способ задаёт лишь их конечное количество, например, экспериментальные данные или результаты сложных вычислений. В связи с этим встаёт вопрос получения промежуточных значений таблично заданной функции. Для их нахождения используют два способа: интерполяция (экстраполяция)*^* Экстраполяция отличается от интерполяции тем, что этим методом пытаются вычислить следующее за последним значение или стоящее до первого заданного, а не промежуточное как при интерполяции. и аппроксимация. При использовании интерполяции считают табличные значения точными и интерполяционный многочлен должен точно проходить через заданные точки (Рис. 11). Аппроксимирующая же функция лишь приближается как можно ближе к заданным точкам, но не обязательно через них проходит (Рис. 13). Это связано с тем, что табличные данные считаются приближенными, и точное прохождение функции через эти точки будет только увеличивать ошибку.

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

Итак, как было сказано выше, задачей интерполяции является поиск такого многочлена, график которого проходит через заданные точки.

Пусть функция y=f(x) задана с помощью таблицы (табл. 1).

Таблица 1

x	x0	x1	x2	…	xn
y	y0	y1	y2	…	yn

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

Необходимо получить многочлен Pn(x) такой, чтобы выполнялось условие:

Pn(xi)=yi. (42)

Для этого зададимся конкретным видом многочлена. Пусть Pn(x) имеет следующий вид:

Pn(xi)=a0+a1x+a2x2+…+anxn. (43)

Для того, чтобы определить коэффициенты a0, a1,… необходимо решить систему из n уравнений с n неизвестными:

(44)

Полином с коэффициентами, полученными путём решения системы (44) называют интерполяционным полиномом Лагранжа и обозначают Ln(x). Решение системы (44) весьма трудоёмко, поэтому интерполяционный полином Лагранжа представляют в виде линейной комбинации многочленов степени n:

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

. (45)

Необходимо, чтобы каждый многочлен li(x) обращался в нуль во всех узлах интерполяции, за исключением i-го, в котором он равен 1 (Рис. 12). Если li(x) удовлетворяет таким условиям, то в i-ом узле интерполяции многочлен Ln(x) примет значение yi, что удовлетворяет условию поставленной задачи. Таким условиям удовлетворяет многочлен вида:

. (46)

Таким образом, интерполяционный многочлен Лагранжа можно представить следующей общей формулой:

(47)

Аппроксимация

Размещено на http://www.allbest.ru/

Размещено на http://www.allbest.ru/

В задачу аппроксимации входит нахождение такой функции y=f(x), что расстояния между заданными точками yi и значениями f(xi) были минимальными (Рис. 13). Обозначим отклонение:

еi=yi-f(xi)

В качестве оценки общего отклонения кривой f(x) от табличных данных (Табл. 1) можно было бы взять сумму отклонений еi, но отклонения могут быть разными по знаку и, не смотря на большие еi их сумма может быть близка к нулю. Очевидно, что необходимо брать сумму абсолютных значений отклонений, но на практике неудобно пользоваться этой функцией, поэтому в качестве критерия оценки отклонения кривой берут сумму квадратов отклонений:

. (48)

Для определения функции f(x) необходимо, во-первых, задать её общий вид, например, f(x)=ax+b, во-вторых, подставив f(x) в (48) и минимизировав у, найти коэффициенты (a и b). Такой метод определения коэффициентов для функции f(x) называется методом наименьших квадратов. Наиболее часто встречающиеся виды функции f(x) для метода наименьших квадратов приведены в таблице 2. Формула y=f(x) называется эмпирической формулой или уравнением регрессии y на x.

Таблица 2

Общий вид функции	Аналитическая формула	Вид регрессии
y=f(x,a,b)	y=ax+b	линейная
y=f(x,a,m)	y=axm	степенная
y=f(x,a,m)	y=aemx	показательная
y=f(x,a,b)		дробно-линейная
y=f(x,a,b)		логарифмическая
y=f(x,a,b)		гиперболическая
y=f(x,a,b)		дробно-рациональная

Рассмотрим подробнее метод наименьших квадратов на примере линейной регрессии, т.е. общий вид функции такой: f(x)=ax+b. Требуется найти методом наименьших квадратов коэффициенты a и b. Для определения минимального у необходимо приравнять нулю частные производные этой функции по параметрам a и b. Для случая линейной регрессии формула (48) приводится к следующему виду:

. (49)

Условие экстремума запишется так:

. (50)

Отсюда получаем систему из двух уравнений с двумя неизвестными.

(51)

(52)



(53)

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

Степенная регрессия

Прологарифмируем степенную зависимость при условии, что x>0 и a>0.

В общем виде линейную зависимость будем записывать так:

Y=AX+B. (54)

Для линейной регрессии Y=y A=a, X=x, B=b. Для регрессии степенной:

.

То есть для того, чтобы воспользоваться формулами (53) необходимо вместо подставлять , вместо - . Тогда по первой формуле определяющей значение параметра A получим значение m, а по второй - значение . Следовательно, остаётся небольшое преобразование для получения a: .

Перепишем формулы (53) в общем виде с учётом (54).

. (55)

Показательная регрессия

Сделаем следующее преобразование:

.

Тогда (a>0, y>0).

Дробно-линейная регрессия.

(y?0).

Логарифмическая регрессия

(x>0).

Гиперболическая регрессия

(x?0).

Дробно-рациональная регрессия

(x?0, y?0).

Численные методы решения дифференциальных уравнений

Постановка задачи

Дано дифференциальное уравнение и начальное условие, то есть поставлена задача Коши:

(2.1.1)

Требуется отыскать интегральную кривую, удовлетворяющую поставленной задаче Коши с помощью метода Рунге-Кутты четвертого порядка с автоматическим выбором шага на отрезке . Задачу можно решить аналитически, найдя решение дифференциального уравнения и подставив в него начальное условие, тем самым, отыскав требуемую интегральную кривую. Но для нас интерес представляет решение данной задачи с применением численного метода, а конкретнее - метода Рунге-Кутты 4-го порядка с автоматическим выбором шага, то есть численное решение. Автоматический выбор шага - необходимое условие адекватного поведения программы при резко изменяющихся функциях, задающих интегральную кривую, позволяющее отразить все моменты в поведении интегральной кривой и добиться высокой точности.

Метод Эйлера

Метод Эйлера для решения начальной задачи (2.1.1) был описан Эйлером в 1768 году. Этот метод весьма прост. Его глобальная погрешность имеет вид , где - постоянная, зависящая от задачи, и - максимальная длина шага. Если желательно, скажем, получить 6 точных десятичных знаков, то требуется, следовательно, порядка миллиона шагов, что не слишком удовлетворительно. С другой стороны, еще со времен Ньютона известно, что можно найти гораздо более точные методы, если не зависит от , то есть если мы имеем задачу (2.1.1), решаемую квадратурой

. (2.2.1)

В качестве примера можно рассмотреть первую квадратурную формулу Гаусса, также называемую «правилом средней точки»:

(2.2.2)

где и - граничные точки подинтервалов, на которые разбит интервал интегрирования. Известно, что оценка глобальной погрешности этой формулы имеет вид . Таким образом, если желаемая точность составляет 6 десятичных знаков, ее обычно можно получить приблизительно за 1000 шагов, то есть этот метод в тысячу раз быстрее. Поэтому Рунге поставил следующий вопрос: нельзя ли распространить этот метод на исходную задачу Коши? Первый шаг длины должен иметь вид

. (2.2.3)

Но какое значение взять для ? За неимение лучшего естественно использовать один малый шаг метода Эйлера длины . Тогда из предыдущей формулы получим:

(2.2.4)

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

(2.2.5)

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

(2.2.6)

Вычитая из последнего равенства предыдущее, получим для погрешности первого шага выражение

(2.2.7)

Таким образом, если все частные производные второго порядка ограничены, то

.

Чтобы получить приближенное значение решения исходной задачи в конечной точке , будем применять формулы (2.2.4) последовательно к интервалам . Приведенные выше формулы являются усовершенствованным методом Эйлера. Для вычислений с высокой точностью, однако, следует пользоваться другими методами, одним из которых как раз является метод Рунге-Кутты.

Общая формулировка методов Рунге-Кутты

Рунге и Хойн построили новые методы, включив в указанные формулы один или два добавочных шага по Эйлеру. Но именно Кутта сформулировал общую схему того, что теперь называется методом Рунге-Кутты.

Пусть - целое положительное число (число стадий, этапов) и - вещественные коэффициенты. Тогда метод

 (2.3.1)

называется -стадийным явным методом Рунге-Кутты для исходной задачи Коши (2.1.1)

Обычно коэффициенты удовлетворяют условиям

. (2.3.2)

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

Метод Рунге-Кутты имеет порядок , если для достаточно гладких задач (2.1.1) справедливо неравенство

, (2.3.3)

то есть ряды Тейлора для точного решения и для совпадают до члена включительно.

После статьи Бутчера вошло в обычай символически представлять метод (2.3.1) по средствам следующей таблицы:

Обсуждение методов порядка 4

Подойдем теперь вплотную к определению 4-стадийных методов Рунге-Кутты (2.3.1) с таким расчетом, чтобы они имели порядок 4. Для этого необходимо вычислить производные порядков 1, 2, 3 и 4 от при и сравнить их с производными точного решения. Теоретически при известных правилах дифференциального исчисления это совершенно тривиальная задача. Однако с использованием (2.3.2) получаются следующие условия:

Эти вычисления очень утомительны и емки. Их громоздкость очень быстро растет для более высоких порядков.

Лемма 1.

Если

(2.4.2)

то уравнения d), g) и h) являются следствием остальных.

Доказательство.

Покажем это для g). C помощью уравнений c) и e) получим:

Для уравнений d) и h) процедура аналогична.

Покажем, что в нашем случае условие

является и необходимым.

Лемма 2.

При (2.4.2) следует из уравнений (2.4.1) и уравнений (2.3.2).

Для доказательства потребуется следующая лемма 3.

Лемма 3.

Пусть и суть 3x3-матрицы, такие что

, (2.4.3)

тогда либо , либо , где .

Доказательство.

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

Докажем теперь предыдущую лемму. Введем величины для . Итак, надо доказать, что . Введем теперь матрицы

(2.4.4)

Перемножение этих матриц с использованием условий (2.4.1) дает

, (2.4.5)

причем

Далее последний столбец не может быть нулевым, так как из того, что , следует

в силу условия h). Таким образом, из последней леммы следует, что . Последнее тождество вытекает из равенства , которое является следствием условий a) и b).

Теорема.

Если выполнены предположения , то уравнения (2.4.1) эквивалентны следующим:

(2.4.6)

Доказательство.

Из j) и h) следует, что

. (2.4.7)

Отсюда, в частности, вытекает, что в силу k) .

Решение уравнений (2.4.6). Уравнения a)-e) и k) выражают тот факт, что коэффициенты и являются весами и узлами квадратурной формулы четвертого порядка при и . В силу (2.4.7) возможны следующие четыре случая:

1) . (2.4.8)

Тогда уравнения a)-e) образуют невырожденную линейную систему для определения . Эта система имеет решение:

Остальные три случая с двойными узлами основаны на правиле Симпсона:

2) ;

3) ;

4) .

После того, как выбраны и , получаем из уравнения j), и тогда два уравнения f) и i) образуют линейную систему для определения и .

Определитель этой системы

,

согласно (2.4.7) не равен нулю. Наконец, из того, что находим , и .

Особенно популярными стали два варианта, которые выбрал Кутта в 1901 году. Это случай 3) при и случай 1) при . Оба метода обобщают классические квадратурные формулы, сохраняя их порядок. Первый из них более популярен, однако второй более точен.

Правило 3/8

Классический метод Рунге-Кутты

«Оптимальные» формулы

Предпринималось много исследования, чтобы выбрать возможно «лучшие» из множества различных формул Рунге-Кутты 4-го порядка.

Первой попыткой в этом направлении был очень популярный метод, который в 1951 году предложил Гилл. Он преследовал цель уменьшить на сколько возможно количество требуемой машинной памяти («регистров»). Этот метод широко использовался на первых компьютерах в пятидесятых годах и представляет поэтому исторический интерес. Гилл заметил, что больше всего машинной памяти нужно при вычислении , когда «требуются регистры для хранения в какой-либо форме» величин

.

Ясно, что для третьей стадии будет достаточно трех регистров, если подлежащие хранению величины линейно зависимы, то есть если

уравнение интеграл касательный секущий

.

Гилл заметил, что это условие удовлетворяется для методов типа 3), если . Получающийся метод можно в таком случае переформулировать следующим образом:



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

Условия порядков для методов Рунге-Кутты

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

Так как явные методы Рунге-Кутты являются частным случаем неявных, то можем выписать условия, при которых метод имеет заданный порядок.

Метод

(где на свободных местах должны стоять нули) имеет порядок , если удовлетворяется уравнение

(2.6.1)

для каждого дерева с корнем и не более чем с разветвлениями Дж. Холл, Дж. Уатт «Современные численные методы решения обыкновенных дифференциальных уравнений», М., Мир, 1979, стр. 77..

При эти условия, обеспечивающие порядок 4, и соответствующие деревья имеют следующий вид:

	(2.6.2)
	(2.6.3)
	(2.6.4)
	(2.6.5)
	(2.6.6)
	(2.6.7)
	(2.6.8)
	(2.6.9)

Заметим, что для меньших значений мы берем соответствующее подмножество этих условий, а для меньших оставляем лишь некоторые из указанных членов.

Из (2.9) видим, что действительно необходимо 4 этапа, так как если бы их было меньше, то был бы опущен единственный член в левой части этого уравнения. Для явных методов в общем случае выполняется неравенство . Фактически (для тех значений, для которых это известно) минимальное значение для данного указано в следующей таблице:

	1	2	3	4	5	6	7
	1	2	3	4	6	7	9

Общие классы методов с этими значениями и легко найти в случае .

Для :

0
	1

Это известный метод Эйлера.

Для :

Это однопараметрическое семейство имеет требуемый порядок для любого ненулевого значения .

Для имеется три семейства, из которых первые два таковы:

Каждое из них имеет один параметр . Третье семейство имеет в качестве параметров и , причем

.

Вывод методов с более сложен, но его можно упростить, положив

(2.6.10)

(что влечет равенство ), так как это позволяет опустить уравнения (2.6.3), (2.6.5), (2.6.8) и (2.6.9). Интересно также, что (2.6.10) является следствием (2.6.2) - (2.6.9).

План вывода конкретного метода этого порядка можно выполнить при условии, что не возникает несовместных систем.

Шаг 1. Выбираем значения , и полагаем .

Шаг 2. Из (2.6.2), (2.6.3), (2.6.4) и (2.6.6) находим .

Шаг 3. Из уравнения (это уравнение есть разность уравнений (2.6.5) и (2.6.7)) находим .

Шаг 4. Из (2.6.10) находим .

Шаг 5. Вычисляем .

В случае шаг 2 приводит к выбору и при условии, что , . В частности, имеем известный метод:

Оценка погрешности и сходимость методов Рунге-Кутты

Со времен работы Лагранжа и особенно Коши всякий установленный численно результат принято сопровождать надежной оценкой погрешности. Лагранж дал известные оценки погрешности многочленов Тейлора, а Коши вывел оценки для погрешности метода ломаных Эйлера. Через несколько лет после первых успехов методов Рунге-Кутты также пришел к заключению, что для этих методов нужны оценки погрешностей «Между тем еще нет доказательства, что эти приближенные методы сходятся, или, что практически важнее, нет критерия, определяющего, сколь малым надо сделать шаги, чтобы достичь предписанной точности» - так писал Рунге в 1905 году..

Строгие оценки погрешности

Способ, которым Рунге получил оценку погрешности, делаемой на одном шаге («локальной погрешности»), может быть описан следующим образом. Для метода порядка рассмотрим локальную погрешность

(2.7.1)

и воспользуемся ее тейлоровским разложением:

, (2.7.2)

где и . Явное вычисление дает выражение вида

, (2.7.3)

где и содержат частные производные до порядков и соответственно. Далее поскольку , имеем . Таким образом, если ограничены все частные производные до порядка включительно, имеем и . Следовательно, существует постоянная такая, что и

. (2.7.4)

Бибербах использовал несколько иной подход. Запишем

(2.7.5)

и воспользуемся тейлоровскими разложениями

(2.7.6)

Для векторных функций эти формулы справедливы покомпонентно (возможно, с различным ). В силу условий порядка первые члены разложения (2.6.5) по степеням обращаются в нуль. Таким образом, справедлива следующая теорема.

Теорема.

Если метод Рунге-Кутты (2.3.1) имеет порядок и если все частные производные до порядка включительно существуют и непрерывны, то локальная погрешность метода (2.3.1) допускает следующую строгую оценку:

, (2.7.7)

или

. (2.7.8)

Продемонстрируем этот результат, применяя к скалярному дифференциальному уравнению первый метод Рунге-Кутты (2.2.4), который имеет порядок . Дифференцируя (2.1.1), получим

. (2.7.9)

Вторая производная величины имеет вид

Если условия теоремы выполнены, то легко видеть, что выражения (2.7.9) и (2.7.10) ограничены постоянной, которая не зависит от , что и дает оценку (2.7.8).

Главный член погрешности

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

Теорема.

Если метод Рунге-Кутты имеет порядок и если непрерывно дифференцируема раз, то для главного члена погрешности имеем:

. (2.7.11)

(2.7.12)

Оценка глобальной погрешности

Глобальной (накопленной) погрешностью Хайрер Э., Нёрсетт С., Ваннер Г. «Решение обыкновенных дифференциальных уравнений. Нежесткие задачи», М., Мир, 1990, стр. 169. называется погрешность численного решения после выполнения нескольких шагов. Пусть мы имеем некоторый одношаговый метод, с помощью которого при заданных начальных данных и длине шага мы определяем численное решение , аппроксимирующее . Воспользуемся обозначениями Хенричи для этого процесса:

, (2.7.13)

и назовем функцией приращения для данного метода.

Оценивание глобальной погрешности методами a) и b)

Тогда численное решение в точке получается с помощью пошаговой процедуры

и наша задача состоит в оценке глобальной погрешности

(2.7.15)

Эта оценка находится простым способом: локальные погрешности переносятся в конечную точку и затем складываются. Этот «перенос погрешностей» можно выполнить двумя разными способами:

a) перенося погрешность вдоль кривых точных решений; этот способ может дать хорошие результаты, если известны хорошие оценки распространения погрешности для точных решений.

b) перенося погрешность -го шага посредством выполнения шагов численного метода; этот способ использовали в своих доказательствах Коши (1824) и Рунге (1905), он легко обобщается на многошаговые методы.

В обоих случаях оценим сначала локальные погрешности:

. (2.7.16)

Займемся теперь оценкой перенесенных погрешностей .

a) Теорема.

Обозначим окрестность точки , где - точное решение уравнения

.

Пусть в справедливы оценки локальных погрешностей (2.7.16) и выполнено одно из условий:

или . (2.7.17)

Тогда имеет место следующая оценка глобальной погрешности (2.7.15):

, (2.7.18)

где ,

и достаточно мало для того, чтобы численное решение оставалось в .

Доказательство.

При оценка (2.7.18) переходит в .

. (2.7.19)

Подставляя в неравенство

выражение (2.7.18) с учетом (2.7.16) и принимая во внимание, что , приходим к такому неравенству:

.

Выражение в квадратных скобках мажорируется следующими интегралами:

, (2.7.20)

. (2.7.21)

Отсюда вытекает справедливость оценки (2.7.18).

b) При втором способе переноса погрешностей рассмотрим кроме (2.7.14) еще одно численное решение, значения которого в соседних узлах связаны равенством

.

Оценим норму разности через . Для формулы метода Рунге-Кутты запишем в следующих обозначениях:

Вычитая из этих формул соответствующие формулы (2.3.1), получим для норм разностей такие оценки:

Оценивание римановых сумм методом a) и b)

Пусть - постоянная Липшица для функции и пусть . Тогда функция приращения для метода (2.3.1) удовлетворяет неравенству

, (2.7.22)

где

. (2.7.23)

Из (2.7.22) получаем искомую оценку:

, (2.7.24)

и с её помощью оценку перенесенных погрешностей вместо оценки (2.7.19).

Предположим, что для начальных значений, лежащих на точном решении, локальная погрешность удовлетворяет оценке

(2.7.25)

и что в окрестности решения функция приращения удовлетворяет неравенству

. (2.7.26)

Тогда для глобальной погрешности (2.7.15) справедлива следующая оценка:

, (2.7.27)

где .

Оптимальный выбор шага

Предположим, что при интегрировании от точки до точки с шагом погрешность приближенно равна . Так как это соответствует росту погрешности со скоростью, приблизительно равной , то , где - функция, определяющая шаг. Положим и получим оценку интеграла, который приближенно равен полной погрешности:

С другой стороны, затраты будут пропорциональны числу шагов, которое приближенно равно



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

В современных программах Амоносов А.А., Дубинский Ю.А., Копченова Н.В. «Вычислительные методы для инженеров», М., Высшая школа, 1994, стр. 445., реализующих методы Рунге-Кутты, обязательно используется некоторый алгоритм автоматического изменения шага интегрирования. Интуитивно ясно, что на участках плавного изменения решения счет можно вести с достаточно крупным шагом. В то же время, на тех участках, где происходят резкие изменения поведения решения, необходимо выбирать мелкий шаг интегрирования. Обычно начальное значение шага задет пользователь или оно определено в программе. Далее шаг интегрирования изменяется в соответствии с величиной, получаемой в ходе вычисления оценки локальной погрешности.

Существует достаточно много способов оценки локальной погрешности, среди которых так называемое правило Рунге. Однако в моей программе я использовал самый простой и в то же время эффективный способ оценки локальной погрешности, который описан в разделе 3.1. «Описание программы Ilya RK-4 версия 1.43». Этот метод базируется на удвоении или делении пополам длины шага в зависимости от отношения локальной погрешности и максимально локальной допустимой погрешности .

Список использованной литературы

1. Калиткин Н.Н. Численные методы. - М.: Наука, 1978. - 512с.

2. Самарский А.А. Введение в численные методы. - М.: Наука, 1982. - 272с.

3. Баранов А.В., Рябчук Э.В. Численные методы в инженерных задачах. - Волгоград.: Политехник, 1988. - 128с.

4. Амоносов А.А., Дубинский Ю.А., Копченова Н.В. «Вычислительные методы для инженеров», М., Высшая школа, 1994, 544с.

5. Хайрер Э., Нёрсетт С., Ваннер Г. «Решение обыкновенных дифференциальных уравнений. Нежесткие задачи», М., Мир, 1990, 512с.

6. Холл Д., Уатт Д. «Современные численные методы решения обыкновенных дифференциальных уравнений», М., Мир, 1979, 312с.

Размещено на Allbest.ru