close

Вход

Забыли?

вход по аккаунту

?

FreeBASIC8

код для вставкиСкачать
Пособие для начинающих программировать на языке высокого уровня FreeBASIC. Работа состоит из нескольких глав. В главе 7 "Решение нелинейных уравнений" – продолжение пособия, которое будет интересно для учащихся школ, студентов институтов, а также пр
 1
Евгений Рыжов, инженер
Программирование на языке FreeBASIC
пособие для начинающих
Фрагмент 8. Решение нелинейных уравнений.
Пособие для начинающих программировать на языке высокого уровня FreeBASIC.
Работа состоит из нескольких глав. В главе 8
"
Решение нелинейных уравнений
" –
продолжение пособия, которое будет интересно для учащихся школ, студентов институтов, а также преподавателей.
Мостовая пусть качнётся, как очнётся,
Пусть начнётся, что ещё не началось,
Вы рисуйте, вы рисуйте, вам зачтётся,
Что гадать нам удалось, не удалось.
"Живописцы", Окуджава Булат,
песня из фильма "Покровские ворота".
2
Не то, чтобы очень сильный стимул для "творчества", но, тем не менее, оправдание процесса –
зачтется!
Булат Шалвович Окуджава
(1924 —
1997) —
советск
ий и российский поэт, композитор, литератор, прозаик и сценарист. При рождении назван родителями Дорианом, в честь Дориана Грея. Автор около двухсот авторских и эстрадных песен, написанных на собственные стихи, один из наиболее ярких представителей жанра а
вторской песни в 1950
—
1980
-
е годы.
"Покровские ворота"
-
советский двухсерийный художественный фильм, снятый Михаилом Козаковым в 1982 году. Литературной основой фильма послужила одноименная пьеса Леонида Зорина. В фильме прозвучали три песни Булата Окудж
авы, фрагменты фильма с которыми можно посмотреть по следующим ссылкам:
"Живописцы"
-
Живописцы окуните ваши кисти...
http://pesnifilm.ru/load/pok
rovskie_vorota/zhivopiscy_quot_pokrovskie_vorota_quot/61
-
1
-
0
-
259
"Песенка об Арбате"
-
Ах, Арбат, мой Арбат, ты мое призвание...
http
://
pesnifilm
.
ru
/
l
oad
/
pokrovskie
_
vorota
/
arbat
_
quot
_
pokrovskie
_
vorota
_
quot
/61
-
1
-
0
-
257
"Часовые любви"
-
Часовые любви на Смоленской стоят...
http
://
pesnifilm
.
r
u
/
load
/
pokrovskie
_
vorota
/
chasovye
_
ljubvi
_
quot
_
pokrovskie
_
vorota
_
quot
/61
-
1
-
0
-
258
На мой взгляд, "Покровские ворота" очень милый и добрый фильм!
1.
Почтовый ящик
Попытаюсь очень коротко ответить на два, пришедшие за последний месяц, письма:
1.1.
Пись
мо от читателя NiV:
Уважаемый NiV спрашивает: нет ли ошибки
в программе "TEST" от 10.10.2012 для определения машинных констант?
Программа приведена в брошюре "Фрагмент 5. Некоторые штришки".
Напомню, при выполнении программа "TEST" на экран дисплея выводя
тся значения "машинного эпсилон" (EPS -
точность численного типа данных):
для типа чисел Single (32 бита) -
5.960464e
-
008
для типа чисел Double (64 бита) -
1.110223e
-
016
Если принять определение: "машинным эпсилон"
называется минимальное положительное чис
ло, которое при сложении с числом 1.0 дает результат, больший, чем 1.0. Очевидно, что исходя из приведенного определения, полученные программой 3
значения, должны быть увеличены в два раза, т.к. они соответствуют максимальному положительному числу, которое п
ри сложении с числом 1.0 дает результат равный 1.0. Замечание следует принять к сведению: "машинное эпсилон" —
наиболее важная характеристика вещественной арифметики. Так в программах из Библиотеки численного анализа НИВЦ МГУ (и не только этой библиотеки) принята константа:
SYS05
1
\
53
-
максимальное число, представимое на данной машине...
Полезную информацию можно найти на сайте:
http://www.softelectro.ru/ieee754.html
IEEE 754 -
стандарт двоичной арифм
етики с плавающей точкой.
Полное название стандарта в ассоциации IEEE:
IEEE Standard for Binary Floating
-
Point Arithmetic
IEEE Стандарт для двоичной арифметики с плавающей точкой
Название стандарта в комиссии IEC:
IEC 60559:1989 Binary floating
-
point arith
metic for microprocessor systems
IEC 60559:1989 Двоичная арифметика с плавающей точкой для микропроцессорных систем
Справка:
Институт инженеров по электротехнике и электронике (
IEEE
-
Institute of Electrical and Electronics Engineers) —
международная нек
оммерческая ассоциация специалистов в области техники, мировой лидер в области разработки стандартов по радиоэлектронике и электротехнике.
Международная электротехническая комиссия (
IEC
-
International Electrotechnical Commission) —
международная некоммер
ческая организация по стандартизации в области электрических, электронных и смежных технологий. Некоторые стандарты МЭК разрабатывает совместно с Международной организацией по стандартизации (ISO).
1.2.
Письмо от читателя AlS:
Уважаемый AlS спрашивает:
не слишком ли "много лирики" содержит брошюра?
Имеется в виду брошюра "Фрагмент 7. Аппроксимация. Сплайны".
Действительно, брошюра кроме листингов программ содержит "внесюжетные рассуждения", но очень сложно указать объем лирических отступлений, который переводит "отчет о научно
-
технической работе" в категорию "мемуары"...
http://www.tdocs.su/13075
Отчет о научно
-
исследовательской работе (НИР) по ГОСТ 7.32
-
2001
Так вот, будем считать, что брошюры включают матери
алы, связанные с выполненной НИР, которые по многим причинам не могли быть включены в основную часть. Но они не собраны в отдельное приложение, а разнесены по всему тексту. Полагаю, это простительное нарушение стандарта. Будем считать, что "отступление от темы" являются своеобразным "педагогическим приемом"!
При помощи авторских отступлений писатели и поэты (избави Бог –
не претендую на эти звания!) как бы приоткрывают завесу над своими мыслями и чувствами, заставляя читателей 4
задуматься о таких непреходящи
х ценностях, как любовь к родине, к людям, уважение, доброта, смелость и самопожертвование...
В этой брошюре предлагаю, например, следующее "авторское отступление"
-
картину русского художника Н.П.Богданова
-
Бельского, написанную им в 1895 году под названи
ем "Устный счёт. В народной школе С.А.Рачинского".
"Устный счёт. В народной школе С.А.Рачинского"
На картине изображена деревенская школа XIX века во время урока арифметики при решении дроби в уме. Учитель —
реальный человек, Сергей Александрович Рачи
нский, ботаник и математик, профессор Московского университета. На волне народничества в 1872 Рачинский вернулся в родное село Татево, где создал школу с 5
общежитием для крестьянских детей, разработал уникальную методику обучения устному счету, прививая дер
евенским ребятишкам его навыки и основы математического мышления. Эпизоду из жизни школы с творческой атмосферой, царившей на уроках, и посвятил своё произведение Богданов
-
Бельский, сам в прошлом ученик Рачинского.
На картине Богданова
-
Бельского изображен
напряженный момент решения в уме дроби:
(10^2 + 11^2 + 12^2 + 13^2 + 14^2)/365
Попробовал и я... Признаюсь –
результат более чем печальный! Пробуйте и вы.
Один из вариантов вычислений, в котором предлагается упростить числитель выражения, по
-
иному сгруппи
ровав его слагаемые, такой:
10^2 + 11^2 + 12^2 + 13^2 + 14^2 = 10^2 + 12^2 + 14^2 + 11^2 + 13^2 =
4*(5^2 + 6^2 + 7^2) + 11^2 + (11 + 2)^2 =
4*(25 + 36 + 49) + 121 + 121 + 44 + 4 =
4*110 + 242 + 48 = 440 + 290 = 730.
Почти уверен, что читатели смогут еще б
олее "упростить" задачу!
Ведь такое стремление присуще всем настоящим инженерам. Напомню, что само слово "инженер" впервые появилось, если верить Анатолию Александровичу Вассерману
(журналист, политический консультант, многократный победитель интеллектуаль
ных телеигр), в 1170
-
ом году, когда английский город Дарем нанял себе в качестве архитектора некоего " Ricardus ingeniator'а", то есть "Ричарда выдумывателя". А в русский язык слово пришло из Франции. Французский вариант слова в XV веке также возник на баз
е латинского ingeneum —
"остроумная выдумка" (в более позднем варианте французское слово "инженьё" переводится как "ухищряться", "проявлять изобретательность").
Монумент Группа инженеров
6
Посмотрите на людей, представленных на монументе "Мемориала при
нца Альберта" (Albert Memorial, монумент в Кенсингстонском парке Лондона). Завидуете? Я –
тоже!
Памятник был открыт королевой Викторией в честь своего мужа Альберта, скончавшегося в 1861 от тифа. Монумент был спроектирован Георгом Гилбертом Скоттом в неого
тическом стиле. Памятник открыт в 1875 и обошелся в 120000 фунтов стерлингов.
Сегодня в нашей стране слово "инженер" означает просто человек с техническим образованием, создатель "материальных средств" для достижения поставленных целей, разработчик способ
ов (технологий) изготовления этих средств
. Как
-
то очень банально звучит, хотя без особого труда можно заметить, непрерывно нарастает разрыв между гуманитарной и технической культурами
(один из факторов разрушая государства). Причем технари традиционно прис
лушиваются к деятелям искусства куда внимательнее, нежели те —
к инженерным и научным творцам. Более того, "инженер или физик, не знакомый с сонетами Шекспира или офортами Гойи, встречается несравненно реже, чем артист или писатель, не только не владеющий школьной математикой, но и считающий правилом хорошего тона сомнение в шарообразности Земли".
2.
Вместо введения. Два математика и один поэт.
Алгебра (от араб. аль
-
джабр —
восполнение) —
раздел математики, который можно грубо охарактеризовать как обобщ
ение и расширение арифметики.
Слово "алгебра" также употребляется в названиях различных алгебраических систем. В более широком смысле под алгеброй понимают раздел математики, посвящённый изучению операций над элементами множества объектов произвольной прир
оды, обобщающий обычные операции сложения и умножения чисел, т.е. выполняемых не только над числами, но и над другими математическими объектами.
2.1.
Математик Нильс Генрих Абель (1802
-
1829)
В Королевском парке в Осло стоит скульптура сказочного юноши,
попирающего двух поверженных чудовищ
: по цоколю идет надпись "ABEL".
Что же символизируют чудовища? Первое из них, несомненно –
алгебраические уравнения 5
-
й степени. Еще в последних классах школы Абелю показалось, что он нашел формулу для их решения, под
обную тем, которые существуют для уравнений степени, не превышающей четырех. Никто в провинциальной Норвегии не смог проверить доказательство. Абель сам нашел у себя ошибку, он уже знал, что не существует выражения для корней в радикалах. Тогда Абель не зн
ал, что итальянский математик П. Руффини опубликовал доказательство этого утверждения, содержащее, однако, пробелы.
Нильс Генрих Абель
7
К тому времени Абель был уже студентом университета в Осло (тогда Христиании). Он был совершенно лишен средств к сущес
твованию, и первое время стипендию ему выплачивали профессора из собственных средств. Затем он получил государственную стипендию, которая позволила ему провести два года за границей. В Норвегии были люди, которые понимали, сколь одарен Абель, но не было та
ких, кто мог бы понять его работы. Будучи в Германии. Абель так и не решился посетить К. Гаусса.
2.2.
Математик Эварист Галуа (1811
-
1832)
Этот человек прожил всего двадцать лет
, всего пять лет из них занимался математикой. Математические работы, обессм
ертившие его имя, занимают чуть более 60 страниц
. В 15 лет Галуа открыл для себя математику и с тех пор, по словам одного из преподавателей, "был одержим демоном математики". Юноша отличался страстностью, неукротимым темпераментом, что постоянно приводило его к конфликтам с окружающими, да и с самим собой. Галуа не задержался на элементарной математике и мгновенно оказался на уровне современной науки. Ему было 17 лет, когда его учитель Ришар констатировал: "Галуа работает только в высших областях математики
". Ему было неполных 18 лет, когда была опубликована его первая работа. И в те же годы Галуа два раза подряд не удается сдать экзамены в Политехническую школу, самое престижное учебное заведение того времени. В 1830 г. он был принят в привилегированную Выс
шую нормальную школу, готовившую преподавателей. За год учебы в этой школе Галуа написал несколько работ; одна из них, посвященная теории чисел, представляла исключительный интерес.
Одной из самых важных задач теории алгебраических уравнений в XVII
-
XVIII
вв. было отыскание формулы для решения уравнения 5
-
й степени. После бесплодных поисков многих поколений алгебраистов усилиями французского ученого XVIII века Ж. Лагранжа (1736
-
1813), итальянского ученого П. Руффини (1765
-
1822) и норвежского математика Н. Абеля в конце XVIII –
начале XIX веков было доказано, что не существует формулы, с помощью которой можно выразить корни любого уравнения 5
-
й степени через коэффициенты уравнения, используя лишь арифметические операции и извлечение корней. Эти исследования были завершены работами Э. Галуа, теория которого позволяет для любого уравнения определить, выражаются ли его корни в радикалах.
Эварист Галуа
О великом французском математике можно прочитать в статье:
Соловьев Ю.П. Эварист Галуа.
http://kvant.mccme.ru/1986/12/evarist_galua.htm
Научно
-
популярный физико
-
математический журнал "Квант" #12
\
1986
8
2.3.
Поэт Владимир Семёнович Высоцкий (1938 —
1980)
25 января 2013 года широко отмечалось 75
-
л
етие Владимира Высоцкого
—
русского советского актера, поэта и автора и исполнителя песен, автора прозаических произведений, лауреата Государственной премии СССР (1987, посмертно). Его знает каждый в нашей стране (думаю и многие в мире). Хотите узнать, как
ое стихотворение написал восьмиклассник Володя 8 марта 1953? Докладываю –
первое стихотворение Высоцкого: "Моя клятва"
(на смерть И.В. Сталина}:
Владимир Высоцкий
Опоясана трауром лент,
Погрузилась в молчанье Москва,
Глубока ее скорбь о вожде,
Сердце бо
лью сжимает тоска...
...
В эти скорбно
-
тяжёлые дни
Поклянусь у могилы твоей
Не щадить молодых своих сил
Для великой Отчизны моей.
И, как мне кажется, Владимир Семёнович исполнил свою клятву!
Очевидно, что события в стране 1972 –
1974 годов не могли никог
о оставить равнодушными –
они породили смутное ощущение, что сломалась какая
-
то очень важная деталь в государственной машине!
Нефтяной кризис 1973 стал своего рода прародителем нынешнего положения России, так как сложившаяся в мире ситуация способствовала увеличению объема поставок нефти из СССР на Запад. Позднее эту "традицию" переняла и расширила Россия. Именно тогда и было положено начало зависимости России от нефтедолларов и "нефтяной трубы".
Послушайте песни Высоцкого:
http://www.youtube.com/watch?v=jrVQaMf5AdQ
Высоцкий: "Чужая колея" (1972)
http://www.youtube.com/watch?v=8Mc5vaZKSX0
Высоцкий: "Погоня" (1974)
http://www.youtube.com/watch?v=mH8_ZISWvLc
Высоцкий: "Чудной дом" (1974)
3.
Решение нелинейных уравнений и систем нелинейных уравнений
Нет спора -
существующие стандартные программные средства существенно облегчают
решение математических задач, но иллюзорна возможность обойтись без знания основ методов вычислений и вообще математики. Уже отмечалось, что имеется обширная учебная литература по методам численного анализа, причем ее "звездные часы" приходились на 40
-
е -
60
-
е годы, когда вычислительная техника находилась на стадии "взросления"
, считалось "хорошим тоном" приводить доступные описания 9
идей методов. При этом строго сочеталась математическая культура изложения с доступностью для заинтересованного читателя, при
водились примеры ручного счета и обсуждались проблемы, которые могут возникнуть при машинном счете (точность, скорость сходимости, затраты машинного времени, памяти). В последующие годы учебная литература резко дифференцировалась: либо издания для "узких с
пециалистов", либо "поваренная книга".
Ну, это так, к слову...
Напомним азы.
Алгебраическое уравнение, уравнение, получающееся при приравнивании двух алгебраических выражений. Алгебраическое выражение, выражение, составленное из букв и цифр, соединённых з
наками действий сложения, вычитания, умножения, деления, возведения в целую степень и извлечения корня (показатели степени и корня должны быть постоянными числами). Алгебраическое выражение называется рациональным относительно некоторых букв, в него входящ
их, если оно не содержит их под знаком извлечения корня. Трансцендентные функции, аналитические функции, не являющиеся алгебраическими. Простейшими примерами Трансцендентные функции служат показательная функция, тригонометрические функции, логарифмическая функция.
К уравнениям, для которых известны аналитические решения, относятся алгебраические уравнения, не выше четвертой степени: линейное уравнение, квадратное уравнение, кубическое уравнение и уравнение четвертой степени. Алгебраические уравнения высших
степеней в общем случае аналитического решения не имеют, хотя некоторые из них можно свести к уравнениям низших степеней. В общем случае, когда аналитического решения найти не удается, применяют численные методы.
Численные методы не дают точного решения, а только позволяют сузить интервал, в котором лежит корень, до определенного заранее заданного значения.
3.1.
Решение нелинейных уравнений
В этом разделе рассмотрим некоторые методы решения уравнений имеющих вид:
F(x) = 0
.
Первая проблема, с которой пр
иходится сталкиваться, это проблема "отделения корней", т.е. нахождение некоторого отрезка [A, B], который содержит корень Z уравнения. Нахождение такого (таких) отрезков может быть задачей более сложной, чем вычисление корня.
3.1.1.
Нахождение отрезка с
одержащего корень
В качестве примера рассмотрим подпрограмму из Библиотеки численного анализа НИВЦ МГУ с сайта: http://num
-
anal.srcc.msu.ru/lib_na/libnal.htm
и далее:
http://num
-
anal.srcc.msu.ru/lib_na/cat/zf/zf18r.htm
Подпрограмма: ZF18R
-
выделение инте
рвала, на котором вещественная функция меняет знак. Подпрограмма определяет границы интервала (A,
B), на котором вещественная функция y = F(x) меняет знак. На входе в подпрограмму задаются предположительные границы интервала (A,
B), длина которого увеличив
ается в FACT раз до тех пор, пока не будет получен интервал (A,
B), на котором F(x) меняет знак, либо пока не будет превзойдено заданное количество итераций.
10
В программе принято:
F -
имя вещественной подпрограммы
-
функции вычисления F(x);
A, B -
вещественн
ые переменные, содержащие на входе в подпрограмму предположительные границы интервала, а на выходе -
найденные границы интервала, на котором F(x) меняет знак;
FACT -
множитель, определяющий, во сколько раз должен увеличиваться текущий интервал (A,
B) на ка
ждом шаге итерации, FACT
>
1;
N -
заданное максимальное количество увеличений первоначального интервала;
IERR –
код ошибки, обнаруженной в ходе работы подпрограммы:
IERR=65 -
когда границы первоначального интервала совпадают;
IERR=66 -
когда превышен лимит
итераций.
Вызывается внешняя подпрограмма: UTZF10 выдачи диагностических сообщений.
Текст программы на языке FreeBASIC и пример использования приведены ниже.
' P R O G R A M "ZeroF18R"
' 12.08.1994
'
-------------------------------------------------------------
' Выделение интервала изменения знака функции
'
-------------------------------------------------------------
' Ординарная переработка подпрограммы ZF18R из:
' Библиотека численного анали
за НИВЦ МГУ
#Lang "FB" ' режим FreeBASIC
-
совместимости
' Иллюстрирующая
подпрограмма
-
функция
Function F(X As Single) As Single
' X -
аргумент функции;
' F -
значение функции.
F = Sin(X)
End Function
' Подпрограмма вывода диагностических со
общений
Sub UTZF10(IERR As Integer, N As Integer)
' Вывод фатальных и не фатальных ошибок выполнения (заглушка)
' IERR –
код обнаруженной ошибки:
' 65 -
границы первоначального интервала совпадают;
' 66 -
интервал изменения знака не найден за N итераций.
' Print IERR; N
End Sub
' Подпрограмма выделение интервала изменения знака функции
Sub ZF18R(ByRef A As Single, ByRef B As Single,_
FACT As Single, ByRef N As Integer, IERR As Integer)
Dim As Integer I, NREAL Dim As Single F1
, F2, ZERO
ZERO = 0.0: IERR = 0
If A = B Then IERR = 65: GoTo L40
F1 = F(A): F2 = F(B)
For I = 1 To N
NREAL = I
If F1*F2 < ZERO Then GoTo L50
If Abs(F1) < Abs(F2) Then GoTo L20
B = B + FACT*(B -
A)
F2 = F(B)
GoTo L30
11
L20: A = A + FACT*(A -
B)
F1 = F(A)
L30: Next I
IERR = 66
L40: UTZF10(IERR, 18)
L50: N = NREAL
End Sub
' 19.01.2013
'
-------------------------------------------------------------
' Основной модуль тестовой программы
'
-------------------------------------------------------------
'
' F -
имя подпрограммы
-
функции вычисления F(x);
' A, B -
переменные, содержащие границы интервала;
' FACT –
множитель увеличения текущего интерв
ала (A,
B);
' N -
максимальное количество увеличений интервала (A,
B).
'
Dim As Single FACT, A, B
Dim As Integer N, IERR
FACT = 1.1 ' возможны варианты...
N = 50 ' максимальное количество увеличений интервала A = 1.0 ' вероятная нижняя граница интервала
B = 1.1 ' вероятная верхняя граница интервала
Open "CON" For Output As #2
'Open "ZeroF18R.res" For Output As #2
ZF18R(A, B, FACT, N, IERR) ' вызов
подпрограммы
Print #2,
Print #2, Using " A = ###.######## B = ###.########"
; A; B
Print #2, Using " NREAL = ## IERR = ##"; N; IERR
Print #2,
Close #2
Sleep
' Результаты вычисления:
' A = -
0.84481049 B = 1.10000002
' NREAL = 5 IERR = 0
'
3.1.2.
Нахождение нуля вещественной функции методом Брента
В качестве одного из примеров рассмотрим подпрограмму из Библиотеки численного анализа НИВЦ МГУ с сайта: http://num
-
anal.srcc.msu.ru/lib_na/libnal.htm
и далее:
http://num
-
anal.srcc.msu.ru/lib_na/cat/zf/zf10r.htm
Подпрограмма: ZF10R
-
вычисление нуля вещественной функ
ции, меняющей знак на заданном интервале, методом Брента. Подпрограмма вычисляет нуль вещественной функции y
=
F(x), меняющей знак на заданном интервале (A,
B), т.е. F(A)*F(B)<0, используя метод Брента.
В программе принято:
F -
имя вещественной подпрограм
мы
-
функции вычисления F(x) в любой точке интервала (A,
B);
12
A, B -
заданная нижняя и верхняя границы интервала, на котором F(x) меняет знак;
EPS -
первый критерий сходимости: заданная абсолютная погрешность;
NDIG -
второй критерий сходимости: заданное число
значащих цифр вычисления нуля;
ITMAX -
целая переменная, значение которой перед началом работы подпрограммы равно максимальному числу итераций, требуемых для обеспечения сходимости; ее значение в результате работы подпрограммы полагается равным действител
ьному числу итераций, потребовавшихся для обеспечения сходимости в соответствии с заданными критериями;
ROOT -
вещественная переменная, значение которой в результате работы подпрограммы полагается равным вычисленному нулю функции F(x);
IERR –
код ошибки в ходе работы подпрограммы:
IERR=65 -
когда нуль не может быть посчитан за заданного числа итераций;
IERR=66 -
когда функция не меняет знака на заданном интервале.
Используется внешняя подпрограмма UTZF10 выдачи диагностических сообщений.
Замечания по исполь
зованию: При обращении к подпрограмме может быть задан только первый критерий сходимости (тогда NDIG задается равным 0), либо только второй критерий (тогда EPS задается равным 0), либо оба критерия одновременно.
Текст программы на языке FreeBASIC и пример
использования приведены ниже.
' P R O G R A M "ZeroF10R"
' 12.08.1994
'
-------------------------------------------------------------
' Вычисление нуля вещественной функции методом Брента
'
--------
-----------------------------------------------------
' Ординарная переработка программы
#Lang "FB" ' режим FreeBASIC
-
совместимости
' Иллюстрирующая
подпрограмма
-
функция
Function F(X As Single) As Single
' X -
аргумент функции;
' F -
значени
е функции.
F = X*X + 2.0*X -
6.0
End Function
' Подпрограмма вывода диагностических сообщений
Sub UTZF10(IERR As Integer, N As Integer)
' Вывод фатальных и не фатальных ошибок при выполнении (заглушка)
' IERR –
код обнаруженной ошибки:
' 65 -
г
раницы первоначального интервала совпадают;
' 66 -
интервал изменения знака не найден за N итераций.
Print IERR; N
End Sub
'
' Подпрограмма вычисление нуля вещественной функции
Sub ZF10R(A As Single, B As Single, EPS As Single,_
NDIG As In
teger, ByRef ITMAX As Integer,_
ByRef ROOT As Single, ByRef IERR As Integer)
Dim As Integer IO, NSIG, NDMAX
Dim As Single A1, B1, C, D, E, S, P, R, RM, RONE
Dim As Single Q, FA, FB, FC, FBC
Dim As Single ZERO,HALF,ONE,T,TEN, TOL, TEMP, T
HREE, SYS051
ZERO = 0.0: HALF = 0.5: ONE = 1.0: THREE = 3.0: TEN = 10.0
13
SYS051 = 1.192093E
-
07 ' машинное эпсилон
NDMAX = 6: A1 = A: B1 = B: NSIG = NDIG: IO = 0: IERR = 0
'
If NDIG = 0 Then NSIG = NDMAX
T = TEN^(
-
NSIG): FA = F(A): FB = F(B)
If FA*FB < 0 Then GoTo L1
If FA*FB = 0 Then GoTo L19
If FA
*FB > 0 Then GoTo L20
L1: C = A: FC = FA: D = B
-
A: E = D
L2: If ABS(FC) >= Abs(FB) Then GoTo L3
A = B: B = C: C = A: FA = FB: FB = FC: FC = FA
L3: TOL = T
If T < (SYS051+SYS051)*ABS(B) Th
en TOL = (SYS051+SYS051)*ABS(B)
If ABS(FB) <
= Abs(EPS) Then GoTo L16
RM = (C
-
B)*HALF
If ABS(RM) >= TOL Then GoTo L4
GoTo L16
L4: If ABS(E) < TOL Or Abs(FA) <= Abs(FB) Then GoTo L12
S=FB/FA
If A < C Then GoTo L6
If A = C Then GoTo L5
If A > C Then GoTo L6
L5: P = (RM+RM)*S
Q = ONE
-
S
GoTo L7
L6: Q=FA/FC
R=FB/FC
RONE=R
-
ONE
P=S*((RM+RM)*Q*(Q
-
R)
-
(B
-
A)*RONE)
Q=(Q
-
ONE)*RONE*(S
-
ONE)
L7: If P < ZERO Then GoTo L8
If P = ZERO Then GoTo L8
If P > ZERO Then GoTo L9
L8: P = -
P
GoTo L10
L9: Q = -
Q
L10: S = E: E = D
IF(P+P) >= THREE*RM*Q
-
ABS(TOL*Q) And P >= Abs(HALF*S*Q) Then GoTo L11
D = P/Q
GoTo L13
L11: E = RM: D=E
GoTo L13
L12: E = RM: D = E
L13: A = B: FA = FB
If ABS(D) >
TOL Then GoTo L14
TEMP=
-
TOL
If RM > ZERO Then TEMP=TOL
GoTo L15
L14: TEMP = D
L15: IO = IO+1
14
If IO <= ITMAX Then GoTo L18
IERR = 65
L16: ITMAX = IO
If IERR <> 0 Then GoTo L21
L17: ROOT = B: B = B1: A = A1
GoTo
L22
L18: B = B+TEMP: FB = F(B): FBC = FB*FC
If FBC < 0 Then GoTo L2
If FBC = 0 Then GoTo L2
If FBC > 0 Then GoTo L1
L19: If FA = ZERO Then B = A
GoTo L16
L20: IERR=66
IO = 1
GoTo L16
L21: UTZF10(IERR,10)
GoTo L17
L22: End Sub
' 19.01.2013
'
-------------------------------------------------------------
' Основной модуль тестовой программы
'
-------------------------------------------------------------
'
' F -
имя подпрограммы
-
функции вычисления F(x);
' A, B -
переменные, содержащие границы интервала;
' EPS -
абсолютная погрешность вычисления нуля функции;
' NDIG -
заданное число значащих цифр вычисления нуля;
' ITMAX -
максимальное число итераций;
' ROOT -
вычисленный нуль функции F(X);
' IERR –
код обнаруженной ошибки:
' 65 –
исчерпан лимит числа итераций;
' 66 -
функция не меняет знак на заданном интервале.
'
Dim As Single A, B, EPS, ROOT
Dim As Integer NDIG, ITMAX, IERR
A = -
3.3: B = 4.6
EPS
= 0.00001
NDIG = 6
ITMAX = 20
ZF10R(A, B, EPS, NDIG, ITMAX, ROOT, IERR)
Print Using " ROOT = ###.#########"; ROOT
Print Using " ITMAX = ### IERR = ##"; ITMAX; IERR
Sleep
' Результат вычисления:
' ROOT = 1.645750999
' ITMAX = 7 IERR = 0
'
15
3.
1.3.
Нахождение нуля вещественной функции по Форсайту
Как вы уже догадались, речь идет о любимой книге:
Форсайт Дж., Малькольм М., Моулер К., Машинные методы математических вычислений. М., Мир, 1980, 280 с., илл.
Глава 7. Решение нелинейных уравнений
За
дача состоит в том, чтобы найти один или более нулей нелинейной функции т.е. решений уравнения f(х) = 0. Нахождение формул для нулей полиномов было одним из важнейших разделов итальянской математики эпохи Ренессанса. Для полиномов 2
-
й, 3
-
й и 4
-
й степеней е
ще несколько столетий назад были найдены алгоритмы, выражающие корни посредством конечного числа квадратичных или кубичных радикалов и рациональных операций. Но только в тридцатых годах прошлого века Галуа доказал невозможность подобных алгоритмов для поли
номов 5
-
й или более высокой степени, даже если допустить в формулах радикалы с показателем n.
Пункт 7.2. Подпрограмма ZEROIN
Приводимая ниже программа иллюстрирует использование подпрограмму ZEROIN() для простого примера f(x) = x^3 –
2*x -
5. Из прикидоч
ного графика функции f(х) легко видеть, что имеется только один действительный нуль и что он лежит между 2 и 3. Результат: Z = 2.0945514815.
В этом примере мы имеем дело с полиномом, к которому можно было бы применить методы 7.4, однако здесь он использова
н по причинам исторического интереса...
' P R O G R A M "Zeroin_Fors"
' 19.01.2013
'
-------------------------------------------------------------
' Нахождение действительного нуля функции ZEROIN
'
-
------------------------------------------------------------
' Ординарная переработка программ из книги:
' Форсайт Дж., Малькольм М., Моулер К., Машинные методы
' математических вычислений. М., Мир, 1980, 280 с., илл.
#Lang "FB" ' режим FreeBASIC
-
со
вместимости
' Иллюстрирующая
подпрограмма
-
функция
Function F(X As Single) As Single
' X -
аргумент функции;
' F -
значение функции.
F = X*(X*X -
2) -
5
End Function
'
' Подпрограмма
-
функция нахождения нуля
Function ZEROIN (AX As Single, BX As
Single, TOL As Single) As Single
' Hуль функции F(X) вычисляется в интервале (AX, BX)
' Входные данные:
' AX -
левая граница интервала
' BX -
правая граница интервала
' F -
значение вычисленное подпрограммой
-
функцией
' TOL -
ширина и
нтервала неопределенности результата
16
' Выходные данные:
' ZEROIN -
абсцисса близкая к нулю функции
' без проверки предполагается, что F(AX) и F(BX) имеют
' противоположные знаки.
' ZEROIN вычисляет "нуль" в интервале (AX, BX)
' в преде
лах допуска на ошибку 4*EPS*ABS(X) + TOL,
' где EPS -
относительная машинная точность.
' Подпрограмма
-
функция представляет собой слегка
' измененную AЛГOЛ
-
60 процедуру ZERO из книги
' RICHARD BRENT, ALGORITHMS FOR MINIMIZATION WITHOUT
' DERIVATIVES,PRENTICE HALL,INC.(1973).
'
Dim As Single A, B, C, D, E, EPS, FA, FB, FC, TOL1, XM, P, Q, R, S
Dim As Single Sign
' Вычисление относительной машинной точности EPS
EPS = 1
Do
EPS = EPS/2
TOL1 = 1 + EPS
Loop While TOL1 > 1
'
' Присвоение начал
ьных значений интервала
A = AX: B=BX
FA = F(A): FB=F(B)
'
' Начать шаг вычислений
L20: C = A: FC = FA: D = B
-
A: E = D
L30: If Abs(FC) >= Abs(FB)Then GoTo L40
A = B: B = C: C = A
FA = FB: FB = FC: FC = FA
'
' Проверка сходимости процесса
L40:
TOL1 = 2.0*EPS*ABS(B)+0.5*TOL
XM = 0.5*(C
-
B)
If ABS(XM) <= TOL1 Then GoTo L90
'
' Оценка необходимости бисекции
If FB = 0 Then GoTo L90
If Abs(E) < TOL1 Then GoTo L70
If Abs(FA) <= Abs(FB) Then GoTo L70
'
' Возможна ли квад
ратичная интерполяция
If A <> C Then GoTo L50
'
' Линейная интерполяция
S = FB/FA
P = 2.0*XM*S
Q = 1.0
-
S
GoTo L60
'
' Обратная квадратичная интерполяция
17
L50: Q = FA/FC: R = FB/FC: S = FB/FA
P=S*(2.0*XM*Q*(Q
-
R)
-
(B
-
A)*(
R
-
1))
Q=(Q
-
1)*(R
-
1)*(S
-
1)
'
' Выбрать знаки
L60: If P > 0 Then Q = -
Q
P=ABS(P)
'
' Приемлема ли интерполяция
If 2.0*P >= 3.0*XM*Q
-
ABS(TOL1*Q) Then GoTo L70
If P >= Abs(0.5*E*Q) Then GoTo L70
E = D: D=P/Q
GoTo L80
'
' Выполнение бисекции
L70: D = XM: E = D
'
' Завершение шага
L80: A = B
FA = FB
If ABS(D) > TOL1 Then B = B + D
If XM >= 0 Then Sign = Abs(TOL1)
If XM < 0 Then Sign = -
Abs(TOL1)
If Abs(D) <= TOL1 Then B = B + Sign
FB =
F(B)
If FB*(FC/ABS(FC)) > 0 Then GoTo L20
GoTo L30
'
' Завершение подпрограммы
L90: ZEROIN = B
End Function
'
'
-------------------------------------------------------------
' Основной модуль тестовой программы
'
-----------------------------
--------------------------------
'
Dim As Single A, B, Z, TOL
A = 2: B = 3
TOL = 1.0E
-
10
Z = ZEROIN(A, B, TOL)
'Open "CON" For Output As #2
Open "Zeroin_Fors.res" For Output As #2
Print #2,Using " A = #.# B = #.# Z = ##.#########"; A; B; Z
Close #2
Slee
p
'
' Результат вычислений:
' A = 2.0 B = 3.0 Z = 2.094551563
'
18
3.1.4.
Нахождение нулей полиномов
Вследствие того что полиномы являются весьма специальными функциями, можно понять наличие множества алгоритмов для вычисления их нулей. Некоторые из н
их принадлежат к числу старейших алгоритмов численного анализа. С минувших столетий ведут свое происхождение методы Горнера, Греффе и Бернулли; в вычислительную эпоху созданы методы Рутисхаузера, Лемера, Лина, Бэрстоу (Bairstow), Бэрайсса (Bareiss) и други
е.
' P R O G R A M "PoliRoot"
' 12.10.2012
'
-------------------------------------------------------------
' Нахождение комплексных корней действительного полинома
'
--------------------------------
-----------------------------
'
' Коэффициенты полинома записаны в порядке убывания степеней
' Использован Пример "3а" из Библиотеки алгоритмов 1б
-
50б.
'
#Lang "fb" ' FreeBASIC compatibility (default)
' Описатель
внешней
подпрограммы
:
Declare Sub Pol
Rt(NN As Integer, AA() As Single,_
XX() As Single, YY() As Single, MM As Integer)
' Основная
программа
:
Dim As Single XCOF(12), ROOTR(12), ROOTI(12), X
Dim As Integer M, J, ErrCd
Data 3 ' полином 3
-
ой степени
Data 1.0, 0.0, -
2
.0, -
5.0
Read M
For J = 0 To M: Read XCOF(J): Next J
'
PolRt(M, XCOF(), ROOTR(), ROOTI(), ErrCd)
'
Open "CON" For Output As #2
'Open "PoliRoot.res" For Output As #2
Print #2,
Print #2, " ErrorCode ="; ErrCd
Print #2,
Print #2, " ROOT
-
R ROOT
-
I"
Print #2,
For J = 1 To M
Print #2, Using " ###.########## ###.##########"; ROOTR(J); ROOTI(J)
Next J
Print #2,
For J = 1 To M
If ROOTI(J) = 0 Then X = ROOTR(J) ' вычисленный
корень
Next J
Print #2, Using " FR0(X) = ##.###########^^^^"; X^3 -
2
*X -
5
X = 2.0945514815 ' значение корня из книги
Print #2, Using " FK1(X) = ##.###########^^^^"; X^3 -
2*X -
5
Print #2, Using " FK2(X) = ##.###########^^^^"; X*(X*X
-
2) -
5
19
Close #2
Sleep
End
' Содержимое
файла
"PoliRoot.res":
' ErrorCode = 0
' ROOT
-
R ROOT
-
I
' -
1.0472757816 -
1.1359398365
' -
1.0472757816 1.1359398365
' 2.0945515633 0.0000000000
' FR0(X) = 9.12119573391E
-
07
' FK1(X) = 9.12119573391E
-
07
' FK2(X) = 9.12119560326E
-
07
'
Sub Pol
Rt(M As Integer, XCOF() As Single,_
ROOTR() As Single, ROOTI() As Single, ErrCd As Integer)
' 12.08.1994
'
-------------------------------------------------------------
' Вычисление комплексных кор
ней действительного полинома
'
-------------------------------------------------------------
'
' Получена ординарной переработкой подпрограммы "POLRT"
' из пакета "IMSL" математических подпрограмм на Fortran.
' Метод последовательных приближений
Берстоу
-
Ньютона.
' Внешние подпрограммы НЕ ИСПОЛЬЗУЮТСЯ!
'
' Параметры вызова:
' M -
степень полинома (1 <= M <= 12)
' XCOF -
массив коэффициентов полинома
' ROOTR -
массив действительных частей корней
' ROOTI -
массив мнимых частей корней
' ErrCd -
код ошибки:
' 0 -
нормальное завершение
' 1 -
M меньше 1 или больше 12
' 2 -
нулевой старший коэффициент
' 3 -
исчерпаны возможности
'
Dim COF(M+2) As Single
Dim As Double SUMSQ, X, XT, XT2, XO, UX, Y
, YT, YT2, YO, UY
Dim As Double TEMP, FI, DX, DY, XPR, YPR, V, U, ALPH
Dim As Integer I, II, N, NX, NXX, N2, L, ICT, FIT, ITEMP
If M < 1 Or M > 12 Then
ErrCd = 1 ' M меньше 1 или больше 12
Else
If XCOF(0) = 0 Then
ErrCd = 2 ' нуле
вой старший коэффициент
Else
If M = 1 Then
ROOTR(1) = -
XCOF(1)/XCOF(0): ROOTI(1) = 0
Else
'
ErrCd = 0
N = M: NX = N: NXX = N + 1: N2 = 1
20
FIT = 0
For L = 0 To N: COF(L+1) = XCOF(L): Next L
'
Do
II = 0: XO = 0.00500101: YO = 0.01000101
'
L50: II = II + 1 ' следующее начальное приближение
X = XO
XO = -
10.0 * YO
YO = -
10.0 * X
X = XO: Y = YO
'
L60: ICT = 0 ' начало итераций
L65: UX = 0.0: UY = 0.0
V = 0.0
XT = 1.0: YT = 0.0
U = COF(N+1) ' младший коэффици
ент
If U = 0 Then X = 0: NX = NX -
1: NXX = NXX -
1: GoTo L130
For I = 1 To N
TEMP = COF(N
-
I+1)
XT2 = X * XT -
Y * YT: YT2 = X * YT + Y * XT
U = U + TEMP * XT2: V = V + TEMP * YT2
FI = CDbl(I) ' преобразование в формат Double
UX = UX + FI * XT * TEMP: UY = UY -
FI * YT * TEMP
XT = XT2: YT = YT2
Next I
'
SUMSQ = UX * UX + UY * UY
If SUMSQ = 0 Then
If FIT = 0 Then GoTo L50
X = XPR: Y = YPR: GoTo L120
End IF
DX = (V * UY -
U * UX) / SUMSQ: X = X + DX
DY =
-
(U * UY + V * UX) / SUMSQ: Y = Y + DY
If ABS(DY) + ABS(DX) < 0.00000001 Then GoTo L100
'
ICT = ICT + 1 ' следующий шаг итерации
If ICT < 500 Then GoTo L65
If FIT = 1 Then GoTo L100
If II >= 5 Then ErrCd = 3: Exit Do ' 500 итер
. от
5 нач
.
приближ
.
GoTo L50
'
L100: For L = 1 To NXX
TEMP = XCOF(L
-
1)
XCOF(L
-
1) = COF(L)
COF(L) = TEMP
Next L
ITEMP = N: N = NX: NX = ITEMP
If FIT = 0 Then FIT = 1: XPR = X: YPR = Y: GoTo L60
L120: FIT = 0
If ABS(Y) < 0.000001 * ABS(X) Then GoTo L130
21
ALPH = X + X
SUMSQ = X * X + Y * Y
N = N -
2
GoTo L140
'
L130: Y = 0: SUMSQ = 0: ALPH = X
N = N -
1
L140: COF(2) = COF(2) + ALPH * COF(1)
For L = 2 To N
COF(L+1) = COF(L+1) + ALPH * COF(L) -
SUMSQ * COF(L
-
1)
Next L
Do
ROOTR(N2) = X:
ROOTI(N2) = Y
N2 = N2 + 1
If SUMSQ = 0 Then Exit Do
Y = -
Y ' комплексносопряженные корни
SUMSQ = 0
Loop
'
Loop While N > 0
End If
End If
End If
End Sub
3.2.
Решение систем нелинейных уравнений
Очень распространенной вычислительной задачей является нахождение некоторых или всех решений системы из k нелинейных алгебраических или трансцендентных уравнений с k неизвестными, имеющей вид:
f1(x1, x2, ... xk) = 0
f2(x1, x2, ... xk) = 0
...
fk(x1, x2, ... xk) = 0
Как и в одномерном случае,
основная проблема состоит в том, чтобы подойти достаточно близко к желаемому корню r чтобы могла начаться быстрая сходимость. На практике можно подчас, итерируя с мужеством и оптимизмом, найти корень без особых трудностей.
3.2.1.
Простая программа реше
ния системы нелинейных уравнений
В качестве одного из примеров рассмотрим подпрограмму из популярной среди бейсиковедов
книги: Дьяконов В.П. Справочник по алгоритмам и программам на языке бейсик для персональных ЭВМ. М., 1989. 240 с.
В книге указано, что
решение систем нелинейных уравнений возможно следующими мето
д
ами:
22
метод простых итераций;
метод Зейделя;
метод Ньютона
-
Рафсона.
Если известно достаточно хорошее начальное приближение к решению системы уравнений F(х) = 0, то эффективным методом повышения точности является метод Ньютона. Метод Ньютона наиболее распространенный метод
решения систем нелинейных уравнений. Он обеспечивает более быструю сходимость по сравнению с методом простых итераций. В основе метода Ньютона лежит идея линеаризации нелинейных
уравнений системы с использованием разложения каждого уравнения системы в ряд Тейлора. Ниже приводится текст на языке FreeBASIC простой программы, реализующей метод Ньютона:
' P R O G R A M "ZeroSys"
' 19.01.2013
'
-------------------------------------------------------------
' Решение системы нелинейных уравнений методом Ньютона
'
-------------------------------------------------------------
' Ординарная переработка программы Программа 4.20 из кни
ги:
' Дьяконов В.П. Справочник по алгоритмам и программам
' на языке бейсик для персональных ЭВМ. М., 1989. 240 с.
'
#Lang "FB" ' совместимость с FreeBasic
Sub SysUr(X() As Single, F() As Single, N As Integer)
' N -
количество уравнений;
' X() -
аргументы функций;
' F() -
значения функций.
F(1) = X(1)+3*(Log(X(1))/Log(10))
-
X(2)*X(2)
F(2) = 2*X(1)*X(1)
-
X(1)*X(2)
-
5*X(1)+1
End Sub
' 19.01.2013
'
----------------------------------------
---------------------
' Основной модуль тестовой программы
'
-------------------------------------------------------------
'
Dim As Integer N, M, S, I, J, K, R
Dim As Single X1, H, EPS
Read N, M, EPS
Dim As Single X(N), F(N), A(N,N), B(N)
Data 2, 50, 0.
0000001
Data 3.4, 2.2
S = 0
For I = 1 To N: Read X(I): Next I
Open "CON" For Output As #2
'Open "ZeroSys.res" For Output As #2
L1: SysUr(X(), F(), N)
For I = 1 To N: B(I)= -
F(I): Next I
For J = 1 TO N
X1 = X(J): H = EPS*Abs(X1): X(J) = X1 + H
SysUr(X(), F(
), N)
For I = 1 To N
A(I,J) = (F(I)+B(I))/H
23
Next I
X(J) = X1
Next J
S = S+1: If S = M+1 Then Print #2, " Stop S = M": GoTo St
For I = 1 To N
-
1
For J = I+1 To N
A(J,I) = -
A(J,I)/A(I,I)
For K = I+1 To N: A(J,K) = A(J,K)+A(J,I)*A(I,K): Next K
B(J) = B(J)+A(
J,I)*B(I)
Next J
Next I
F(N) = B(N)/A(N,N)
For I = N
-
1 To 1 Step -
1
H = B(I)
For J = I+1 To N: H = H
-
F(J)*A(I,J): Next J
F(I) = H/A(I,I)
Next I
R = 0
For I = 1 To N
X(I) = X(I) + F(I)
If Abs(F(I)/X(I)) > EPS Then R = 1
Next I
If R = 1 Then GoTo L1
Print #2, " Solving:"
For I = 1 To N: Print #2, Using " ###.#########"; X(I): Next I
Print #2, " Iter ="; S
St: Close #2
Sleep
' Результаты вычислений:
'
' Solving:
' 3.487442732
' 2.261628628
' Iter = 18
3.2.2.
Решение системы нелинейных уравнений мет
одом Ньютона
В отличие от программы, приведенной в разделе 3.2.1. здесь приводится полный вариант программы Решение системы нелинейных уравнений методом Ньютона из Библиотеки численного анализа НИВЦ МГУ:
http://num
-
anal.srcc.msu.ru/lib_na/cat/zf/zf31r.htm
ZF31R -
Решение системы N нелинейных уравнений с N неизвестными методом Ньютона.
В ней используются внешние подпрограммы:
AFG5R
-
треугольное разложение матрицы методом Гаусса;
ASG9R
-
решение системы уравнений методом Гаусса;
каждая из которых допускает собственное использование.
24
3.2.2.1. Внешняя подпрограмма AFG5R
Подпрограмма осуществляет треугольное разложение вещественной матрицы общего вида методом Гаусса с выбором ведущего элемента по столбцу. Пример использования подпрограммы приведен ниже:
' P R O G R A M "MatrAFG5R"
' 10.10.2012
'
-------------------------------------------------------------
' Треугольное разложение матрицы методом Гаусса
'
--------------------------------------------------
-----------
' Библиотека численного анализа НИВЦ МГУ
' http://num
-
anal.srcc.msu.ru/lib_na/cat/af/afg5r.htm
' Треугольное разложение вещественной матрицы общего вида
' методом Гаусса с выбором ведущего элемента по столбцу.
'
Dim As Integer M, N, I, J, IERR
M = 4: N = 4
Dim As Single A(N, N)
Dim As Integer NLEAD(N)
Data 7.9, 8.5, 4.3, 3.2
Data 5.6, -
4.8, 4.2, -
1.4
Data 5.7, 0.8, -
3.2, -
8.9
Data -
7.2, 3.5, 9.3, 3.3
For I = 1 To N
For J = 1 To N: Read A(J,I): Next J
Next I
Open "CON" For Output As #2
'Open "MatrAFG5R.res" For Output As #2
AFG5R(A(), M, N, NLEAD(), IERR)
Print #2,
Print #2, " *** A(I,J):"
For I = 1 To N
For J = 1 To N: Print #2, Using " ###.#####"; A(I,J);: Next J
Print #2,
Next I
Print #2,
Print #2, " *** NLEAD(J):"
For J = 1 To N: Print #2, Using " ###.#####"; NLEAD(J);: Next J
Print #2,
Print #2,
Close #2
Sleep
' Результаты
вычислений
:
' *** A(I,J):
' 8.50000 -
4.80000 0.80000 3.50000
' -
0.92941 10.06118 4.95647 -
10.45294
' -
0.50588 -
0.65879 -
9.40171
2.40526
' -
0.37647 -
0.04046 -
0.73072 12.65817
' *** NLEAD(J):
' 2.00000 2.00000 4.00000 4.00000
25
3.2.2.2. Внешняя подпрограмма ASG9R
Подпрограмма предназначена для решения системы Ax = b или ATx = b методом Гаусса с выбором ведущего элем
ента по столбцу. Пример использования подпрограммы приведен ниже:
' P R O G R A M "MatrASG9R"
' 10.10.2012
'
-------------------------------------------------------------
' Решение системы нелинейны
х уравнений методом Ньютона
'
-------------------------------------------------------------
' Библиотека численного анализа НИВЦ МГУ
' http://num
-
anal.srcc.msu.ru/lib_na/cat/as/asg9r.htm
' Решение системы уравнений методом Гаусса с выбором ведущего
' элемента по столбцу. Используется подпрограмма AFG5R –
' подпрограмма треугольной факторизации матрицы.
'
Dim As Integer M, N, I, J, LTR, L, IERR
M = 4: N = 4: LTR = 0: L = 0
Dim As Single A(N, N), B(N)
Dim As Integer NLEAD(N)
Data 7.9, 8.5, 4.3, 3.2
Data 5.6, -
4.8, 4.2, -
1.4
Data 5.7, 0.8, -
3.2, -
8.9
Data -
7.2, 3.5, 9.3, 3.3
Data 7.0, 7.0, 7.0, 7.0
For I = 1 To N
For J = 1 To N: Read A(J,I): Next J
Next I
For J = 1 To N: Read B(J): Next J
Open "CON" For Output As #2
'Open "MatrAFG5R
.res" For Output As #2
ASG9R(A(), M, N, NLEAD(), B(), LTR, L, IERR)
Print #2,
Print #2, " *** A(I,J):"
For I = 1 To N
For J = 1 To N: Print #2, Using " ###.#####"; A(I,J);: Next J
Print #2,
Next I
Print #2,
Print #2, " *** NLEAD(J):"
For J = 1 To N: Print #2, Using " ###.#####"; NLEAD(J);: Next J
Print #2,
Print #2,
Print #2, " *** B(J):"
For J = 1 To N: Print #2, Using " ###.######"; B(J): Next J
Print #2,
Print #2,
Close #2
26
Sleep
' Результаты
вычислений
:
'
' *** A(I,J):
' 8.50000 -
4.80000 0.80000 3.50000
' -
0.92941 10.06118 4.95647 -
10.45294
' -
0.50588 -
0.65879 -
9.40171 2.40526
' -
0.37647 -
0.04046 -
0.73072 12.65817
'
' *** NLEAD(J):
' 2.00000 2.00000 4.00000 4.00000
'
' *** B(J):
' 1.023054
' 0.273777
' -
0
.462958
' -
0.003275
3.2.2.3. Сборка подпрограмм решения системы
Подпрограмма ZF31R
решения системы N нелинейных уравнений с N неизвестными методом Ньютона достаточно полно описана на сайте Библиотеки численного анализа НИВЦ МГУ. Для общего понимания на
сайте приведена книга:
Н.С.Бахвалов. Численные методы. М., 1979.
Текст сборки подпрограмм решения системы на языке FreeBASIC приведен ниже:
' P R O G R A M "ZeroSF31R"
' 17.01.2013
'
-----------------
--------------------------------------------
' Решение системы нелинейных уравнений методом Ньютона
'
-------------------------------------------------------------
' Ординарная переработка программы
#Lang "FB" ' режим FreeBASIC
-
совместимости
' Иллю
стрирующая
подпрограмма
-
функция
Sub SysUr(X() As Single, ALPH() As Single, BETA() As Single)
BETA(1) = -
(X(1) + EXP(X(1) -
1) + (X(2) + X(3))^2 -
27)
BETA(2) = -
(X(1)*EXP(X(2) -
2) + X(3)^2 -
10)
BETA(3) = -
(X(3) + SIN(X(2) -
2) + X(2)^2 -
7)
A
LPH(1, 1) = 1 + EXP(X(1) -
1)
ALPH(1, 2) = 2*(X(2) + X(3))
ALPH(1, 3) = 2*(X(2) + X(3))
ALPH(2, 1) = EXP(X(2) -
2)
ALPH(2, 2) = X(1)*EXP(X(2) -
2)
ALPH(2, 3) = 2*X(3)
ALPH(3, 1) = 0.0
ALPH(3, 2) = COS(X(2) -
2) + 2*X(2)
ALPH(3, 3) =
1.0
End Sub
'
Sub UTZF10(IERR As Integer, N As Integer)
27
' "Заглушка":
' Вывод фатальных и не фатальных сообщений об ошибках
' IERR –
код обнаруженной ошибки:
' 65 –
исчерпан лимит итераций;
' 66 -
якобиан заданной системы имеет особенность.
Print
IERR; N
End Sub
'
' Подпрограмма вывода диагностических сообщений
Sub UTAFSI(NAM1 As Integer, NAM2 As Integer, NAM3 As Integer, NAM4 As Integer, NAM5 As Integer, IERR As Integer)
' "Заглушка":
' Вывод фатальных и не фатальных сообщений об ошибках
' IERR –
код обнаруженной ошибки:
IF IERR = 0 Then GoTo L9
Print IERR; NAM1; NAM2; NAM3; NAM4; NAM5
L9:
End Sub
'
Sub AFG5R(A() As Single, M As Integer, N As Integer,_
NLEAD() As Integer, IERR AS Integer)
'
---------------------------------------
----------------------
' Треугольное разложение матрицы методом Гаусса
'
-------------------------------------------------------------
' LU
-
разложение матрицы общего вида А методом Гаусса
' с выбором ведущего элемента по столбцу.
' M –
число строк м
атрицы
' N –
размер матрицы
' A –
матрица M*N, на входе -
исходная, на выходе -
элементы
' матрицы U и поддиагональные элементы матрицы INVERSE(L)
' NLEAD –
вектор длины N, на выходе вектор перестановок;
' IERR –
код ошибки.
'
' установлен
ие контроля переполнения
Dim As Integer I, J, JJ, K, KP1, L, NM1
Dim As Integer NAM1, NAM2, NAM3, NAM4, NAM5
Dim As Single T
'
' проверка корректности параметров M и N
If N >= 1 And M >= 1 Then GoTo L20
IERR = 65
UTAFSI(NAM1, NAM2, NAM3, NAM4, NAM5, IERR)
GoTo L200
L20:
IERR = 0
NM1 = N -
1
If NM1 < 1 Then GoTo L190
For K = 1 To NM1
KP1 = K + 1
'
' нахождение индекса ведущего элемента L
28
L = K
T = ABS(A(K,K))
For J = KP1 To N
If ABS(A(J,K)) <= T Then GoTo L30
L = J
T = ABS(A(J,K))
L30:
Next J
NLEAD(K) = L
If A(L,K) = 0 Then GoTo L160
'
' перестановка, если требуется
If L = K Then GoTo L50
T = A(L,K)
A(L,K) = A(K,K)
A(K,K) = T
'
' вычисление множителей
L50: T = -
1./A(K,K): J = N -
K
I = (J Mod 5) + K ' целый остаток от деления
If I = K Then GoTo L70
For J
= KP1 To I
A(J,K) = T*A(J,K)
Next J
L70:
I = I + 1
IF I > N Then GoTo L90
For J = I To N Step 5
A(J,K) = T*A(J,K)
A(J+1,K) = T*A(J+1,K)
A(J+2,K) = T*A(J+2,K)
A(J+3,K) = T*A(J+3,K)
A(J+4,K) = T*A(J+4,K)
L80: Next J
L90:
'
' Гауссово исключение
For J = KP1
To N
T = A(L,J)
If L = K Then GoTo L100
A(L,J) = A(K,J)
A(K,J) = T
L100: JJ = N -
K
I = (JJ Mod 4) + K
If I = K Then GoTo L120
For JJ = KP1 To I
A(JJ,J) = A(JJ,J) + T*A(JJ,K)
Next JJ
L120: I = I + 1
If I > N Then GoTo L140
For JJ = I To N Step 4
A(JJ,J) =
A(JJ,J) + T*A(JJ,K)
29
A(JJ+1,J) = A(JJ+1,J) + T*A(JJ+1,K)
A(JJ+2,J) = A(JJ+2,J) + T*A(JJ+2,K)
A(JJ+3,J) = A(JJ+3,J) + T*A(JJ+3,K)
Next JJ
L140:
Next J
GoTo L170
L160: IERR = -
K
L170:
Next K
L190:
NLEAD(N) = N
If A(N,N) = 0 Then IERR = -
N
IF IERR = 0 Then GoTo L200
UTAFSI(NAM1, NAM2, NAM3, NAM4, NAM5, IERR)
L200:
End Sub
'
Sub ASG9R(A() As Single, M As Integer, N As Integer,_
NLEAD() As Integer, B() As Single, LTR As Integer,_
L As Integer, IERR As Integer)
'
-----------------------
--------------------------------------
' Решение системы уравнений методом Гаусса '
-------------------------------------------------------------
' Система уравнений A*X=B или TRANS(A)*X=B решается методом
' Гаусса с выбором ведущего элемента по столбцу
для матрицы
' A общего вида.
' Используются внешние подпрограммы: AFG5R и UTAFSI.
'
Dim As Integer I, J, JJ, K, KB, KP1, NM1
Dim As Integer NAM1, NAM2, NAM3, NAM4, NAM5
Dim As Single S, T
S = 0.0
'
' проверка корректности параметров M и N
If N >= 1 And
M >= 1 Then GoTo L20
IERR = 65
UTAFSI(NAM1, NAM2, NAM3, NAM4, NAM5, IERR)
GoTo L330
L20: IERR = 0
NM1 = N -
1
If L <> 0 Then GoTo L30
AFG5R(A(), M, N, NLEAD(), IERR)
If IERR <> 0 Then UTAFSI(NAM1, NAM2, NAM3, NAM4, NAM5, IERR)
If IERR <= 0 Then GoTo L30
GoTo L330
L30: If LTR <> 0 Then GoTo L190
'
' Решение системы L*Y=B
If NM1 < 1 Then GoTo L100
For K = 1 To NM1
30
JJ = NLEAD(K)
T = B(JJ)
If JJ = K Then GoTo L40
B(JJ) = B(K)
B(K)
= T
L40: JJ = N -
K
I = (JJ Mod 4) + K
If I = K Then GoTo L60
KP1 = K + 1
For JJ = KP1 To I
B(JJ) = B(JJ) + T*A(JJ,K)
Next JJ
L60: I = I + 1
If I > N Then GoTo L90
For JJ = I To N Step 4
B(JJ) = B(JJ
) + T*A(JJ,K)
B(JJ+1) = B(JJ+1) + T*A(JJ+1,K)
B(JJ+2) = B(JJ+2) + T*A(JJ+2,K)
B(JJ+3) = B(JJ+3) + T*A(JJ+3,K)
L70: Next JJ
L90: Next K
'
' Решение системы U*X=Y
L100: For KB = 1 To N
K = N + 1 -
KB
If A(K,K) <> 0 Then GoTo
L120
If B(K) <> 0 Then GoTo L110
B(K) = 1.0
GoTo L130
L110: IERR = 67
UTAFSI(NAM1, NAM2, NAM3, NAM4, NAM5, IERR)
GoTo L330
L120: B(K) = B(K)/A(K,K)
L130: T = B(K)
J = K -
1
I = (J Mod 4)
If I = 0 Then GoTo L
150
For JJ = 1 To I
B(JJ) = B(JJ) -
T*A(JJ,K)
L140: Next JJ
L150: I = I + 1
If I > J Then GoTo L180
For JJ = I To J Step 4
B(JJ) = B(JJ) -
T*A(JJ,K)
B(JJ+1) = B(JJ+1) -
T*A(JJ+1,K)
B(JJ+2) = B(JJ+2) -
T*A(JJ+2,K)
B(JJ+3) = B(JJ+3) -
T*A(JJ+3,K)
L160: Next JJ
L180: Next KB
GoTo L330
'
' Решение системы TRANS(U)*Y=B
31
L190: For K = 1 To N
T = 0
J = K -
1
I = (J Mod 5)
If I = 0 Then GoTo L210
For JJ = 1 To I
T = T + A(JJ,
K)*B(JJ)
L200: Next JJ
L210: I = I + 1
If I > J Then GoTo L230
For JJ = I To J Step 5
T = T + A(JJ,K)*B(JJ) + A(JJ+1,K)*B(JJ+1) + A(JJ+2,K)_
*B(JJ+2) + A(JJ+3,K)*B(JJ+3) + A(JJ+4,K)*B(JJ+4)
L220: Next JJ
L230: If A(K,K) <> 0 Then G
oTo L250
If B(K) -
T <> 0 Then GoTo L240
B(K) = 1.0
GoTo L260
L240: IERR = 67
UTAFSI(NAM1, NAM2, NAM3, NAM4, NAM5, IERR)
GoTo L330
L250: B(K) = (B(K) -
T)/A(K,K)
L260: Next K
'
' Решение системы TRANS(L)*X=Y
If NM1 < 1
Then GoTo L330
For KB = 1 To NM1
K = N -
KB
JJ = N -
K
I = (JJ Mod 5) + K
If I = K Then GoTo L280
KP1 = K + 1
For JJ = KP1 To I
B(K) = B(K) + A(JJ,K)*B(JJ)
L270: Next JJ
L280: I = I + 1
If I > N Then
GoTo L300
For JJ = I To N Step 5
B(K) = B(K) + A(JJ,K)*B(JJ) + A(JJ+1,K)*B(JJ+1) +_
A(JJ+2,K)*B(JJ+2) + A(JJ+3,K)*B(JJ+3) + A(JJ+4,K)*B(JJ+4)
L290: Next JJ
L300: JJ = NLEAD(K)
If JJ = K Then GoTo L320
T = B(JJ): B(JJ) = B(K):
B(K) = T
L320: Next KB
L330:
End Sub
'
Sub ZF31R(ROOT() As Single, N As Integer, EPSX AS Single,_
EPSF AS Single, ITMAX As Integer, ALPH() As Single,_
BETA() As Single, NLEAD() As Integer, IERR As Integer,_
32
ByRef ITER As Integer)
' 19.01.2013
'
-------------------------------------------------------------
' Решение системы N нелинейных уравнений методом Ньютона
'
-------------------------------------------------------------
'
Dim As I
nteger K, I
Dim As Single ERRF, ZERO, ERRX
' COMMON/ZF31RR/ITER
ZERO = 0: IERR=0
For K = 1 To ITMAX
ITER=K
SysUr(ROOT(), ALPH(), BETA())
ERRF = ZERO
For I = 1 To N
ERRF=ERRF+ABS(BETA(I))
M10: Next I
If ERRF <= Abs(EPSF) Then GoTo M
60
ASG9R(ALPH(), N, N, NLEAD(), BETA(), 0, 0, IERR)
If IERR <> 0 Then GoTo M40
ERRX=ZERO
For I = 1 To N
ERRX=ERRX+ABS(BETA(I))
ROOT(I)=ROOT(I)+BETA(I)
M20: Next I
If ERRX <= Abs(EPSX) Then GoTo M60
M30: Next K
IERR=65
GoTo M50
M40: IERR=66
M50: UTZF10(IERR, 31)
M60: End Sub
'
' 19.01.2013
'
-------------------------------------------------------------
' Основной модуль тестовой программы
'
---------------------------------
----------------------------
' Библиотека численного анализа НИВЦ МГУ
' http://num
-
anal.srcc.msu.ru/lib_na/cat/zf/zf31r.htm
' ZF31R -
Решение системы N нелинейных уравнений с N
' неизвестными методом Ньютона.
' Используется внешние подпрограммы:
' AFG5R –
треугольная факторизация матрицы;
' ASG9R -
решение системы методом Гаусса;
' SysUr() -
вычисление левых частей и производных.
' N -
заданное число уравнений системы;
' EPSX -
первый критерий сходимости;
' EPSF -
второй крит
ерий сходимости;
' ITMAX -
максимальному числу итераций;
'
Dim As Integer N, ITMAX
33
N = 3
Dim As Single ROOT(N), ALPH(N, N), BETA(N)
Dim AS Single IERR, EPSX, EPSF
Dim As Integer J, NLEAD(N), ITER
Data 1.0, 1.0, 1.0 ' начальное
приближение
For J = 1 To N: Read ROOT(J): Next J
EPSX = 0.00001
EPSF = 0.00001
ITMAX = 30
ZF31R(ROOT(), N, EPSX, EPSF, ITMAX, ALPH(), BETA(), NLEAD(), IERR, ITER)
Open "CON" For Output As #2
'Open "ZeroSF31R.res" For Output As #2
Print #2, " *** ROOT(J):"
For J = 1 To N: Print #2, Using " ###.#####"; ROOT(J): Next J
Print #2, Using " ITER = ### IERR = ###"; ITER; IERR
Close #2
Sleep
'
' Результаты вычислений:
' *** ROOT(J):
' 1.00000
' 2.00000
' 3.00000
' ITER = 7 IERR = 0
Сборка получилась, действитель
но, могучая! Простите, я ее только немного адаптировал под цели настоящей брошюры...
Общий недостаток всех рассмотренных выше методов решения систем нелинейных уравнений состоит в локальном характере сходимости, затрудняющем их применение в случаях (доста
точно типичных), когда существуют проблемы с выбором начального приближения, обеспечивающего сходимость итерационной вычислительной процедуры.
B этих случаях можно использовать численные методы оптимизации!
Если все сложится, об этих методах поговорим в следующем выпуске.
На этом позвольте закруглиться.
До новых встреч!
Пишите: eugene_r@mail.ru
Документ
Категория
Информатика
Просмотров
163
Размер файла
620 Кб
Теги
программы, freebasic, алгоритмы
1/--страниц
Пожаловаться на содержимое документа