Говорухин В., Цибулин Б. Компьютер в математическом исследовании (1185927), страница 92
Текст из файла (страница 92)
Для иллюстрации разберем нахождение корней кубического полинома и построение карт начальных значений, из которых итерации сходятся к каждому корню. Метод Ньютона описан во всех руководствах по численному анализу и численным методам, см., например, 1575, К78, С841. Пусть задано скалярное уравнение Г(х) = 0 В этом случае расчет приближений к корню ведется по формуле =012... Л.) У(~5 ) Для начала вычислений необходимо начальное приближение лк итерации продолжаются, пока разность хл+г - лл не станет достаточно малой. Для системы нелинейных уравнений Г по х есть векторы, а производная от Г по х — матрица Якоби.
Библиотека МА6 Замечательная коллекция процедур численного анализа на Фортране поддерживается фирмой ХАС ).Ы. В МАТ1АВ на основе атой библиотеки имеется пакет Тоо)Ьох йАб, куда входит более 250 процедур решения разнообразных задач. Для просмотра содержимого нужно обратиться к справке или выполнить команду Ие1р Нлй Пакет состоит из следующих разделов (в сгюбках будем указывать префикс раздела): гз вычисление нулей полиномов (с02); О нахождение корней трансцендентных уравнений и систем (с06); о вычисление рядов, преобразования Фурье и т.д. (с0б); о квадратуры для одномерных и многомерных интегралов (001); о методы решения обыкновенных дифференциальных уравнений (задачи с начальными данными и краевые задачи) (г)02); Решение нелинейных уравнений 489 О интегрирование уравнений в частных производных (((03); О задачи интерполяции (е01); О аппроксимация кривых и поверхностей (е02); О задачи минимизации (е04); О факторизация матриц ((01); О вычисление собственных значений и векторов (Г02); О решение систем линейных алгебраических уравнений ([04); О процедуры из пакета?.АРАСК (г07); О обработка статистических данных (001); О вычисление корреляций и регрессионный анализ (д02) и (003); О генераторы случайных чисел (005); О методы непараметрической статистики (д08); О анализ временных рядов (д13); О методы сортировки (ш01); О аппроксимация специальных функций (001), (013)-(015), (017)-(021).
Для справки по функциям библиотеки нужно знать имя интересующей нас функ- ции, которое состоит из двух частей: префикса раздела (три буквы) и идентифика- тора (три буквы). Например, для получения справки по функции с05пЬГ нужно набрать » Пе)р с05по( Решение системы нелинейных уравнений можно разыскивать при помощи функ- ции с05пЬ[, которая требует только вычисления значений функций правой части. Кроме того, в пакете Тоо1Ьох МА0 имеется функция с05рЬ(, для использования кото- рой необходимо знать якобиан системы (матрицу первых производных от уравне- ний по переменным), Обращение к функции с05пЬГ имеет вид: [х,т,(та(11 - с05пЬГ(топ,хо,хго),Ч[аа() Обязательными параметрами являются имя функции для вычисления правой ча- сти системы гоп и начальное приближение к решению х0. По умолчанию расчет ведется с точностью адгг(ерз), где ерз есть константа МАТ1.АВ для используемого компьютера.
Точность можно изменить, задав параметр х101; целая переменная ) га(1 по умолчанию равна -1. Для расчета итераций применяется вариант метода Нью- тона с численным определением матрицы Якоби направленными разностями. Выходными параметрами являются найденное решение х, значение функции на решении т и переменная ошибки (га(1. Для определения корней подготовим функцию вычисления правой части в файле с05 йв. Обратим внимание, что имя функции в заголовке может быть произволь- ным (например, топ), а функция опознается по имени файла. Кроме того, в списке параметров должна присутствовать размерность вектора неизвестных х, даже если она не используется при вычислении значений функций: » гипс((оп [глт)а01 гоп(п.хлт)ад) т (х(1)*(х(1)"2-3"х(2)"2)+1.
х(2)*(З»хы)'2.х(2)"2)1: 490 Глава 10, Дополнения и примеры Затем обратимся к функции с05пйб с полным набором выходных параметров: » Гх,гига!))=с05пьг("с05 к.(! 11) х 5.0000е-001 8.6603е-001 -1.2992е-012 3 5959е-014 гба1) = О В результате выполнения команды найдено одно пз решений, нулевое значение переменной ) Га)1 подтверждает, что вычисления завершились нормально. Обращение к функции с05рьт выглядит следующим образом: (х Л,Сас.1га1)) = с05рот(тсп.хО.л(о) л Га1!) Здесь входные параметры те же, что и для функции с05пЬ1. К выходным параметрам добавилась матрица Якоби Оас. Для итераций используется стратегия, сочетающая л(етод Ньютона с методом градиентного спуска. В функции гоп вычисляется правая часть системы и матрица Якоби. Для зтого предусмотрена переменная т1ад, которая при значении 2 вычисляет значения матрицы Якоби, а при прочих значениях — правую часть системы: Гппсг, оп (Г.
С а с, Г) ад) 6бсп(п. х, Г ма с, ) об)ас т) ад) т 1г г)ад-2 Г " Гх(!)*(х(1)"2-3*х(2)"2)+1. х(2)*(3*х(1) 2-х(2) "2)1. е) ае Сас =[3*х(1) "2-3*х(2) 2 -6*х(1)*х(2): б*х(1)*х(2) 3*х(1) "2-3*х(2) 2); еп0 Обращение к функции с05рЬГ с начальным значением хО=Г-1 Ц приводит к тому же результату, что и для функции с05пЬГ. Для вычисления другого решения изменим начальное приближение: » х-с05рвб("с05 Г 3",(-1 11); х' апа = -1.0000е 000 1.5601е-013 Пример функции для решения системы нелинейных уравнений методом Ньютона Приведем пример программирования функции для нахождения корней системы нелинейных уравнений.
Реализуем прямой метод Ньютона, который будет работать как для задачи с введенной матрицей Якоби, так и для случая численного определения матрицы Якоби. Далее представлен текст функции пеысопв вместе с необходимыми пояснениями. В заголовке функции указан порядок следования входных и выходных параметров, коммента ии объясняют их назначение: Решение нелинейных уравнений 491 типот!оп [гО.п] " пеысопв(тип, гО, ивах, бе1, Ь. геви)2) Ж Мвигопв — реюение систеиы нелинейных уравнений т(х)-0 Ж тип — иия функции правой части Ж гд - начальное приблииение Ж ивах - предельное число итераций Ж бе) - точность нетода Ньютона Ж П - юаг для вычисления якобиана численно Ж (отсутствие.
() или [] - и. Якоби дана в функции [ип) Ж гези12 — иня файла для записи итераций Ж Припер вызова: Ж [х.п]-пеитопз("с05 т)3".[1 1],10,!е-б.)е-з,'гези12.2х2') Далее идет анализ входных параметров и задание их значений в том случае, если при вводе они были опущены.
Важным является пятый параметр — величина шага Ь. Если этот параметр отсутствует, на его месте находится пара квадратных скобок или он равен нулю, то этим указывается функции пеытопз, что вычисление матрицы Якоби содержится в функции вычисления уравнений [ип. Если необходима запись результатов решения системы в файл, то в качестве третьего и четвертого параметров могут фигурировать пустые квадратные скобки: 1т пагдтп<2.
01зр("Нет начального приблииения").гетигп епб тт пагд1п<3 ) 1зевргу(ивах). ниах=9: епб Ж нет ивах 1! пагд1п<4 ) 1зеврсу(бе1). бе1-вцг2(ерв);епб Ж нет бе1 1[ пагдтп<5. Ь=зцгт(ерь); епб Ж нет Ь 1т пагд1п<5 ) 1вевртубй) ) Ь-О, )асоь-1; Ь-0: Ж Матрица Якоби дана фориулой е)зе ]асоЬ-0; Ж Матрица Якоби рассчитывается епб Начальное приближение может быть задано вектором или матрицей.
В расчетах будем использовать вектор-столбец, так что вначале производится соответствующее преобразование. После этого следует обращение к функции ]асо, которая помещается в том же файле и предназначена для вычисления правой части и матрицы Якоби; пг-1епдтшгО( )): гО=гО(:); Ж Столбец [тгО,)ас]-)асо(тип.гО,пг,)асов.Ь): Если указано имя файла [шестой параметр), то нужно приготовить переменные, в которые будет заноситься информация об итерациях и достигнутой точности (отклонении от нуля): тт пагдтп-б, г- гд' Я-попа(-тг0); епб Теперь запускается цикл для расчета итераций, где на каждом шаге метода Ньютона анализируется матрица Якоби. Вычисления прекращаются, если матрица Якоби содержит нечисловой параметр, бесконечность или если она вырождена; тог и-М пах бе2)ас - бег(бас): тт -апу(!в[1п12е(бас(:))) ) -1вттп!2е(бетбас).
сдвр("Матрица Якоби плоха"): гесигп е1ве1т ббе2]ас — 0). бе12а"-бас\гад: гО-гд+бе12а: 2а92 . Глава тв. дополнения и примеры е!ае Шар("матрица Якоби необратина"): геьцгп епб ЕгаО Лас]-]асо(гцп,аО.па,басоЬ,Ь); теат-погш((ад) то!ге!-погш(бе1(а)/(попа(гО)тере): тт пагдтп-=б. а(п+1,:)ааО': Я=(Я; теат]: епб тт то)ге! < бе1 $ теат <= бе).
сопчегд - 1; Ьгеай: е)аетт п-пшах. сопчегд = 0; Грю птт( цгНет схоаиностит2.0т итерацийтгтп".ошах): епб епб Выход нз цикла происходит по достижении заданной точности или если будет до- стигнуто максимальное число итераций пп)ах. Результат последней итерации пре- вращается в вектор-строку для использования в выходных параметрах; аОгаО'; Если шестым параметром было указано имя файла для записи результатов, то про- исходит сохранение подготовленной информации о расчете: тт пагдтп > 5. Еп,па] - а!ге(г): Гтбйбореп(геац) Ц "н"); тргтптт(гтб.'Итерации Норна ! Решение тгтп'); тог К"1:п тргтптт(ттб, ГЖЗ.ОГ Ж!2.2Е".
Е, й(Е)); тог ш 1:п7. тргтптт(1тб,'$15.3Е т($11.3Е)',... геа1(а(й,ш)).таад(а(й,ш))). епб Грг! птт(ттб. ''тгтп' ): епб тт -сопчегд тргтптт(гтб,''тгтпНет сходиности$2.0т итерацийтг'тп',ошах); епб тс!оае(т!б); епб Функция ]асс предназначена для вычисления значений функции Ецп и матрицы Якоби. Если в функции Гцп не содержится явных формул для матрицы Якоби, то производится расчет с использованием односторонних разностей: Гцпсттоп Ета.дас]-басс(гцп.г,па.Засао,п) тт )асоЬ, Ета.)ас]-теча!(гцп.г. 1); е)ве Гг-теча)(гцп,а): тог и - 1 па бе иегов(пг.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
