close

Вход

Забыли?

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

?

FreeBASIC12

код для вставкиСкачать
Пособие для начинающих программировать на языке высокого уровня FreeBASIC. Работа состоит из нескольких глав. В главе 12 "Численное интегрирование" – продолжение пособия, которое будет интересно для учащихся школ, студентов институтов, а также преп
 1
Евгений Рыжов, инженер
Программирование на языке FreeBASIC
пособие для начинающих
Фрагмент 12. Численное интегрирование
Пособие для начинающих программировать на языке высокого уровня FreeBASIC.
Работа состоит из нескольких глав. В главе 12 "Чис
ленное интегрирование" –
продолжение пособия, которое будет интересно для учащихся школ, студентов
институтов, а также преподавателей.
Толпой огромною стеснилися в мой ум
Разнообразные, удачные сюжеты,
С завязкой сложною, с анализом души
И с патетичною, з
агадочной развязкой.
2
Я думал в "мировой поэме" их развить,
В большом, посредственном иль в маленьком масштабе.
И уж составил план. И, к миросозерцанью
Высокому свой ум стараясь приучить,
Без задней мысли, я к простому пониманью
Обыденных основ стремился в
сей душой.
Но, верный новому в словесности ученью,
Другим последуя, я навсегда отверг:
И личности протест, и разочарованье,
Теперь дешевое, и модный наш дендизм,
И без основ борьбу, страданья без исхода,
И антипатии болезненной причуды!
А чтоб не впасть в абсурд, изгнал экстравагантность...
Очистив главную творения идею
От ей несвойственных и пошлых положений,
Уж разменявшихся на мелочь в наше время,
Я отстранил и фальшь и даже форсировку
И долго изучал без устали, с упорством
Свое, в изгибах разных, внутре
ннее "Я".
Козьма Прутков "Безвыходное положение"
(г. Аполлону Григорьеву, по поводу статей его
в "Москвитянине" 1850
-
х годов)
Козьма Петрович Прутков
—
литературная маска, под которой в журналах "Современник", "Искра" и других выступали в 50
—
60
-
е годы
XIX века поэты Алексей Толстой (наибольший в количественном исчислении вклад), братья Алексей, Владимир 3
и Александр Жемчужниковы (фактически —
их коллективный псевдоним), а также Пётр Ершов.
1.
Почтовый ящик.
Ждал, ждал позитива и дождался!
Читатель S
S пишет
: Хотите, чтобы Ваши брошюры читали школьники?
"Пишите проще и к Вам потянуться люди"!
Уважаемый SS, большое спасибо за мудрый совет!
Пожалуй, Вы правы –
буду стараться... но в меру своих, как Вы понимаете, не безграничных возможностей... Афоризм: "Будь проще и люди к тебе потянутся" приписывают Козьме Пруткову. Поэтому он и попал на "титульный лист" брошюры. Правда, сегодня чаще используют следующую "трансформацию": "Будь проще, и люди на тебе оттянутся"
... "Но, как же быть с идеями?" –
спросите Вы
. И Вам ответят на вопрос –
вопросом
(как в Одессе)
: "А эти идеи принесут позитивные плоды?"
Посмотрим главный справочник:
http://ru.wikipedia.org/wiki/Творчество
Он предлагает следующее определение:
Творчество
—
процесс деятельности, создающий качественн
о новые материальные и духовные ценности или итог создания объективно нового. Основной критерий, отличающий творчество от изготовления (производства) —
уникальность его результата. Результат творчества невозможно прямо вывести из начальных условий. Никто, кроме, возможно, автора, не может получить в точности такой же результат, если создать для него ту же исходную ситуацию. Таким образом, в процессе творчества автор вкладывает в материал некие несводимые к трудовым операциям или логическому выводу возможнос
ти, выражает в конечном результате какие
-
то аспекты своей личности. Именно этот факт придаёт продуктам творчества дополнительную ценность в сравнении с продуктами производства.
Основная загвоздка данной дефиниции заключается в том, как именно определить к
ачественную новизну
, отличить ее от некачественной? При радикальном подходе сегодня можно, конечно, легко избежать этого вопроса, принудительно именуя любую новую структуру "творчеством". Однако при этом утрачивается специфика данного понятия, оно размывае
тся, становясь простым синонимом слова "новое". Некоторое уточнение сделал Сергей Леонидович Рубинштейн
(1889
—
1960) —
выдающийся российский психолог и философ, член
-
корреспондент Академии наук СССР, автор оригинальной интерпретации психологии как научной д
исциплины в ключе марксистской философии. Он, пожалуй, впервые правильно указал на характерные особенности изобретательского творчества: "Специфика изобретения, отличающая его от других форм творческой интеллектуальной деятельности, заключается в том, что оно должно создать вещь, реальный предмет, механизм или приём, который разрешает определенную проблему. Этим определяется своеобразие творческой работы изобретателя: изобретатель должен ввести что
-
то новое в контекст действительности, в реальное протекание
какой
-
то деятельности".
Можно еще вспомнить наиболее известное описание последовательности стадий (этапов) творческого мышления, которое дал англичанин Грэм Уоллес
в 1926 году. Он выделил четыре стадии творческого мышления:
4
1.
Подготовка —
формулирование за
дачи; попытки её решения.
2.
Инкубация —
временное отвлечение от задачи.
3.
Озарение —
появление интуитивного решения.
4.
Проверка —
испытание и/или реализация решения.
Впрочем, это описание не оригинально и восходит к классическому докладу А. Пуанкаре
1908 года. Т
ак вот, из
-
за необходимости соблюсти пункт "Инкубация"
брошюры кажутся избыточно "рыхлыми" –
простите.
Но ответ читателю SS был бы однобоким, если не вспомнить Курта Фридриха Гёделя
(Kurt Friedrich Godel; 1906 -
1978) —
австрийского логика, математика и философа математики, наиболее известного сформулированной и доказанной им теоремой о неполноте. Как свидетельствует переписка К.Геделя, он был самостоятельным мыслителем
-
теистом, и проявлял особый интерес к исследованию онтологического аргумента бытия Божи
я. Важным аспектом философской мысли Геделя является то, что он был не только теистом, но и персоналистом. Он глубоко ощущал негативный пафос материалистического позитивизма, и рассматривал свою теорему как сокрушительное оружие против него. Понимая под "Г
ёделизацией" cведение любого алгоритма к численному алгоритму
, проникнемся важностью проделываемой работы и продолжим нашу тему, но не забудем и нашу "художественную галерею"
–
в этой брошюре пусть будет картина Мстислава Добужинского...
Добужинский Мстис
лав Валерианович (1875
–
1957)
, русский художник. Родился в Новгороде 14 августа 1875 в семье офицера. Отец Добужинского был родом из Литвы, что сказалось на последующей биографии Добужинского. В 1924 художник принял литовское гражданство и уехал из России, в 1925 работал для Рижского театра, с 1926 по 1929 —
для парижского театра Н.Ф.Балиева "Летучая мышь", с 1929 Добужинский —
ведущий художник литовского Государственного театра. В 1935 с группой каунасского театра уехал в Англию, с 1939 жил в США и скончалс
я в Нью
-
Йорке 20 ноября 1957. В творчестве Добужинского особое место занимает тема современного индустриального города, приобретающего в его произведениях фантастически
-
одушевленный облик. Но, пожалуй, не это главное...
Гражданственность искусства Добужин
ского ярко выразилась во время событий революции 1905
-
1907 годов. Художник выступил с серией резких политических карикатур в журналах "Жупел" и "Адская почта". Особенно выразителен представленный ниже рисунок "Умиротворение"
, где изображена Москва, затопле
нная морем крови. В это время Добужинский с группой деятелей "Мира искусства" принимает участие в организации сатирических журналов и выступает в печати с обращением к художникам, призывая их работать для народа. А я вспоминаю Москву недавнего прошлого
и с
лова сказанные тогда: "Господин президент! Почему Вы стреляли в свой народ? Верховный Совет плохой? Так это Ваш Верховный Совет, это парламент Вашего народа, и другого народа у Вас не будет. Потому что Вас выбрали президентом всея Руси, а не только своей к
оманды"...
5
Умиротворение. 1905. (рисунок для журнала "Жупел", 1905, № 2)
2.
Введение.
Ну, а теперь следуя указаниям читателя SS
"конспективно" о численном интегрировании. Для развитого дифференциального или интегрального исчисления характерно, что после строгого обоснования своих основных понятий при помощи предельного перехода они дают возможность решать разнообразнейшие задачи при помощи простого алгоритма чисто алгебраического характера (в том смысле, что сам этот алгоритм уже не содержит в явном
виде предельных переходов). Благодаря этому современные способы вычисления с дифференциалами и интегралами успешно соединяют в себе строгую логическую обоснованность с простотой и наглядностью.
Нахождение числового значения интеграла является одной из на
иболее распространенных вычислительных процедур при проведении научных исследований и инженерных расчетов. Если заданная подынтегральная функция имеет первообразную в виде сравнительно простого аналитического выражения, то эта процедура не вызывает осложне
ний и состоит в проведении вычислений, связанных с подстановкой числовых значений пределов интегрирования в формулу Ньютона —
Лейбница. Однако в случае сложной первообразной ее использование для вычисления может быть не всегда рационально. Если же интеграл
"неберущийся" или подынтегральная 6
функция задана табличным способом
(например, в виде результатов экспериментальных измерений), то аналитическое выражение интеграла вообще отсутствует. В таких случаях приходится проводить численное интегрирование, под кот
орым понимают процедуру нахождения приближенного значения интеграла методами вычислительной математики.
Интегрирование прослеживается ещё в древнем Египте, примерно в 1800 году до Р.Х., Московский математический папирус демонстрирует знание формулы объёма
усечённой пирамиды. Первым известным методом для расчёта интегралов является метод исчерпывания Евдокса
(примерно 370 до Р.Х.), который пытался найти площади и объёмы, разрывая их на бесконечное множество частей, для которых площадь или объём уже известны
. Этот метод был подхвачен и развит Архимедом, и использовался для расчёта площадей парабол и приближенного расчёта площади круга. Аналогичные методы были разработаны независимо в Китае в 3
-
м веке н.э. Лю Хуэйем, который использовал их для нахождения площа
ди круга. Этот метод впоследствии использовали Цзу Чунчжи и Цзу Гэн для нахождения объёма шара.
Евдокс Книдский (Eudoxus; ок. 408 год до н.э. —
ок. 355 год до н.э.)
—
древнегреческий математик, механик и астроном. О жизни Евдокса известно немного. Родилс
я в Книде, на юго
-
западе Малой Азии. Учился медицине, потом математике (у пифагорейца Архита в Италии), затем присоединился к школе Платона в Афинах. Около года провёл в Египте, изучал астрономию в Гелиополе. Позднее Евдокс переселился в город Кизик на Мра
морном море, основал там собственную математико
-
астрономическую школу, читал лекции по философии, астрономии и метеорологии.
Около 368 года до н.э. Евдокс вместе с частью учеников вернулся в Афины. Умер в родном Книде, окружённый славой и почётом. Диоген Л
аэртский сообщает, что у Евдокса были три дочери и сын по имени Аристагор. Евдокс получил фундаментальные результаты в различных областях математики. Например, при разработке своей астрономической модели он существенно продвинул сферическую геометрию. Одна
ко особенно большое значение имели созданные им две классические теории:
Общая теория отношений
Метод исчерпывания
Следующий крупный шаг в исчисление интегралов был сделан в Ираке, в XI веке, математиком Ибн ал
-
Хайсамом
(известным как Alhazen в Европе), в
своей работе "Об измерении параболического тела" он приходит к уравнению четвёртой степени. Решая эту проблему, он проводит вычисления, равносильные вычислению определённого интеграла, чтобы найти объём параболоида. Используя математическую индукцию, он с
мог обобщить свои результаты для интегралов от многочленов до четвёртой степени. Таким образом, он был близок к поиску общей формулы для интегралов от полиномов, но он не касается любых многочленов выше четвёртой степени.
7
Следующий значительный прогресс в исчислении интегралов появится лишь в XVI веке. В работах Кавальери
с его методом неделимых, а также в работах Ферма
, были заложены основы современного интегрального исчисления. Дальнейшие шаги были сделаны в начале XVII века Барроу и Торричелли
, которые
представили первые намеки на связь между интегрированием и дифференцированием.
Всем известно обозначение операции интегрирования, однако Ньютон использовал в качестве символа интегрирования значок квадрата (перед обозначением функции или вокруг него), но
это обозначение не получило широкого распространения. Современное обозначение неопределённого интеграла было введено Лейбницем в 1675 году. Он образовал интегральный символ из буквы S —
сокращения слова лат. Summa. Современное обозначение определённого ин
теграла, с указанием пределов интегрирования, были впервые предложены Жаном Батистом Жозефом Фурье
в 1819
-
20 годах.
Численное дифференцирование
, т.е. совокупность методов вычисления значения производной дискретно заданной функции в этой брошюре рассматрив
ать не будем –
для нас это уже слишком простая задача! В основе численного дифференцирования лежит аппроксимация функции интерполяционным многочленом, от которой аналитически берется производная...
3.
Знакомство с численными методами интегрирования.
Чи
сленное интегрирование (историческое название: квадратура) —
вычисление значения определённого интеграла (как правило, приближённое). Под численным интегрированием понимают набор численных методов отыскания значения определённого интеграла. Квадратура (лат
. Quadratura, придание квадратной формы) —
математический термин, первоначально обозначавший нахождение площади заданной фигуры или поверхности. В дальнейшем смысл термина постепенно менялся. Задачи квадратуры послужили одним из главных источников возникно
вения в конце XVII века математического анализа. Сразу оговоримся –
будем рассматривать численное интегрирование только для одномерного случая
(функций одной переменной). Полученные результаты при необходимости несложно использовать для нахождения кратных интегралов.
3.1.
Квадратурные формулы Ньютона
-
Котеса.
Для большинства функций поиск первообразной нереален
, не говоря уж о табличных функциях, и потому возникает задача численной оценки интеграла с какой
-
то точностью. Большинство известных квадратурных
формул строится на задании системы точек с абсциссами A + kH (k = 0, 1,..., m), где H=(B -
A)/N. Если пределы интегрирования A и B входят в состав узлов некоторой квадратурной формулы, то ее называют квадратурой открытого типа (в противном случае -
замкну
того типа). Создаваемые таким образом формулы называются формулами Ньютона
-
Котеса. П.Л. Чебышевым
была предложена группа квадратурных формул с неравноотстоящими узлами, но c одинаковым весом слагаемых. Квадратурные формулы Гаусса
отличаются от формул Ньюто
на
-
Котеса и Чебышева тем, что здесь подлежат определению и узлы интегрирования, и весовые коэффициенты.
8
В брошюре "Фрагмент 7. Аппроксимация. Сплайны" при обсуждении главы "Сплайн
-
аппроксимация, интерполяция и экстраполяция" из книги: Дьяконов В.П. Справо
чник по алгоритмам и программам на языке бейсик для персональных ЭВМ. 1989, была упомянута косинус
-
функции Френеля
С(х) (Fresnel cosine integral). Здесь предлагается вычислить значения этой функции как интеграла
. Соответствующая программа на языке FreeBASI
C приведена ниже:
' P R O G R A M "Quadra0"
' 19.01.2013
'
-------------------------------------------------------------
' Интегрирование функции методом прямоугольников
'
--------------------------
-----------------------------------
' на примере вычисления косинус
-
функции Френеля С(х)
' Calculates the Fresnel cosine integral C(x).
' Рассматривается метод прямоугольников (нулевой порядок)
' В программе обозначено:
' X -
переме
нная интегрирования
' A -
нижняя граница интегрирования (фиксирована)
' B -
верхняя граница интегрирования (подвижна)
' N -
количество вычисляемых точек функции
' H -
шаг интегрирования
' C -
вычисленное значение инт
еграла
'
#Lang "FB" ' режим совместимости с FreeBASIC
' Подпрограмма
-
функция вычисления подынтегральной функции
Function Fun(X As Single) As Single
Const As Single Pi = 3.14159265
' X -
аргумент
подынтегральной
функции
Fun = Cos(Pi*X*X/2)
End Func
tion
'
' Подпрограмма интегрирования методом прямоугольников
Sub QRECT(V As Integer, A As Single, B As Single, _
N As Integer, ByRef C As Single)
Dim As Single X, H, S
Dim As Integer I
H =(B
-
A)/N: S = 0 ' Величина шага и накопленной суммы
If V = 1 Then X = A ' Вариант "левые прямоугольники"
If V = 2 Then X = A+H/2 ' Вариант "середины прямоугольников" If V = 3 Then X = A+H ' Вариант "правые прямоугольники"
For I = 1 To N ' Выполнить для N точек
S = S + Fun(X) ' Приба
вить значение функции
X = X +H ' Прибавить шаг к аргументу
Next I ' Следующая точка
C = S*H ' Вычисление "площади"
End Sub
'
' Основной модуль программы
Dim As Integer N, I, J, JJ, V, VV
Dim As Single A, B, HB, C, CT(4), S2
Dim As Single X, Y(128), H, Ymin, Ymax
9
' Табличные значения функции Френеля:
' C(X) X
Data 0.492344225871 ' 0.5
Data 0.779893400377 ' 1.0
Data 0.445261176040 ' 1.5
Data 0.488253406075 ' 2.0
'
' Вывод графика подынтегрально
го выражения Fun(X)
A = 0: B = 2: N = 128 ' Границы интервала и количество точек
H =(B
-
A)/N: X = A
Ymin = 1E38 ' Начальное "минимальное" значение
Ymax = -
1E38 ' Начальное "максимальное" значение
For I = 0 To N ' Перебор всех точек
Y(I) = Fun(X) ' Вычисление значения функции
If Ymin > Y(I) Then Ymin = Y(I) ' Заменить
MIN
If Ymax < Y(I) Then Ymax = Y(I) ' Заменить
MAX
X = X+H
Next I
Screen 12 ' Графический экран 640x480
View (1, 1)
-
(478, 478), , 15 ' Квадратное "физическое" окно
Window (0, Ymin)
-
(N, Ymax) ' "Логическое" окно вывода
Line (0, 0)
-
(128, 0), 7, , &HAAAA ' Линия оси X
For I = 0 To N ' Перебор аргумента
PSet(I, Y(I)), 11 ' Вывод точек графика
Next I ' Следующая точка
' Вв
од табличных значений функции
JJ = 4
For J = 1 To JJ
Read CT(J)
Next J
VV = 3 ' Число вариантов вычисления
N = 8 ' Количество вычисляемых точек
Open "CON" For Output As #2
'Open "Quadra.res" For Output As #2
Print #2,
Print #2, " Int
egration by the Method of Rectangles"
Print #2, Using " Ymin = ##.##### Ymax = ##.#####"; Ymin; Ymax
Print #2,
For V = 1 To 3 ' Вариант вычисления интеграла
If V = 1 Then Print #2, " Option Left
-
rectangles"
If V = 2 Then Print #2, " Option Mid
-
rectangles" If V = 3 Then Print #2, " Option Right
-
rectangles"
A = 0: B = 0.5: HB = 0.5
S2 = 0 ' Сумма квадратов отклонений
For J = 1 To JJ ' Перебор точек
QRECT(V, A, B, N, C) ' Метод прямоугольников
Print #2, Using " X = ##.### C(
X) = ##.#######"; B; C
S2 = S2 + (C
-
CT(J))^2 ' Добавка квадрата отклонения
B = B + HB ' Следующая точка
Next J
S2 = S2/N
Print #2, Using " S2 = ###.#######"; S2
10
Print #2,
Next V
Close #2
Sleep
' Результаты
вычислений
:
' Integration by the M
ethod of Rectangles
' Ymin = -
0.99942 Ymax = 1.00000
' Option Left
-
rectangles
' X = 0.500 C(X) = 0.4945276
' X = 1.000 C(X) = 0.8382922
' X = 1.500 C(X) = 0.6309028
' X = 2.000 C(X) = 0.4886461
' S2 = 0.0047348
' Option Mid
-
rectangles
'
X = 0.500 C(X) = 0.4924418
' X = 1.000 C(X) = 0.7819479
' X = 1.500 C(X) = 0.4426246
' X = 2.000 C(X) = 0.4879028
' S2 = 0.0000014
' Option Right
-
rectangles
' X = 0.500 C(X) = 0.4897701
' X = 1.000 C(X) = 0.7132922
' X = 1.500 C
(X) = 0.2701754
' X = 2.000 C(X) = 0.4886462
' S2 = 0.0043872
Ниже приведена программа для сравнения "точностных" характеристик
трех простых методов интегрирования:
1.
Метод прямоугольников (нулевой порядок)
2.
Метод трапеций (первый порядок)
3.
Метод п
арабол (второй порядок)
' P R O G R A M "Quadra1"
' 20.02.2013
'
-------------------------------------------------------------
' Интегрирование функции одной действительной переменной
'
-------------
------------------------------------------------
' на примере вычисления косинус
-
функции Френеля С(х)
' Calculates the Fresnel cosine integral C(x).
' Рассматриваются методы интегрирования Ньютона
-
Котеса:
' 1. Метод прямоугольников (нуле
вой порядок)
' 2. Метод трапеций (первый порядок)
' 3. Метод парабол (второй порядок)
'
' В программе обозначено:
' X -
переменная интегрирования
' A -
нижняя граница интегрирования (фиксирована)
' B -
верхняя гран
ица интегрирования (подвижна)
' N -
количество вычисляемых точек функции
' H -
шаг интегрирования
' S -
вычисленное значение интеграла
'
11
#Lang "FB" ' режим совместимости с FreeBASIC
' Вычисление
подынтегральной
функции
:
Function Fu
n(X As Single) As Single
Const As Single Pi = 3.14159265
' X -
аргумент
подынтегральной
функции
Fun = Cos(Pi*X*X/2)
End Function
'
' Подпрограмма интегрирования методом прямоугольников:
Sub QRECT(A As Single, B As Single, N As Integer, ByRef C As Single)
Dim As Single X, H, S
Dim As Integer I
H =(B
-
A)/N: X = A+H/2: S = 0 ' Шаг
, начальная
точка
, сумма
For I = 1 To N ' Выполнить для N точек
S = S + Fun(X) ' Прибавить значение функции
X = X + H ' Прибавить шаг к аргуме
нту
Next I ' Следующая точка
C = S*H ' Вычисление "площади"
End Sub
'
' Подпрограмма интегрирования методом трапеций:
Sub QTRAP(A As Single, B As Single, N As Integer, ByRef C As Single)
Dim As Single X, H, S
Dim As Intege
r I
H = (B
-
A)/N: X = A
S = (Fun(A)+Fun(B))/2 ' Крайние
точки
For I = 1 To N
-
1 ' Внутренние
точки
X = X + H
S = S + Fun(X) Next I
C = S*H End Sub
'
' Подпрограмма интегрирования методом парабол:
Sub QPARB(A As Single, B As Single, N As Integer, ByRef C As Single)
Dim As Single X, H, S
Dim As Integer I
H = (B
-
A)/N: X = A
S = (Fun(A)+Fun(B))/2 ' Крайние
точки
I = 0
Do
I = I + 1: X = X + H
S = S + 2*Fun(X)
If I >= N
-
1 Then Exit Do
I = I + 1: X = X + H
S = S + Fun(X)
Loop
C = 2*H*S/3
End Sub
12
'
' Основной модуль программы
Dim As Integer N, J, JJ, M, MM
Dim As Single A, B, HB, C, CT(4), S2
' Табличные значения функции Френеля:
' C(X) X
Data 0.492344225871 ' 0.5
Data 0.779893400377 ' 1.0
Data 0.445261176040 ' 1.5
Data 0.4882534060
75 ' 2.0
' Ввод табличных значений
JJ = 4
For J = 1 To JJ
Read CT(J)
Next J
MM = 3 ' Число методов интегрирования N = 8: A = 0: HB = 0.5 ' Количество точек, нижняя граница, шаг
Open "CON" For Output As #2
'Open "Quadra.res" For Output As #2
For M = 1 To MM
Print #2,
If M = 1 Then
Print #2, " Integration by the Method of Rectangles"
B = 0.5: S2 = 0 ' Верхняя граница и сумма квадратов
For J = 1 To JJ ' Перебор точек
QRECT(A, B, N, C) ' Метод прямоугольников
Print #2, Using " X = ##.### C(X) = ##.#######"; B; C
S2 = S2 + (C
-
CT(J))^2 ' Добавка квадрата отклонения
B = B + HB ' Следующая точка
Next J
S2 = S2/N
Print #2, Using " S2 = ##.########"; S2
Print #2,
EndIf
If M = 2 Then
Print #2, " Integration by the Method of Trapezoids"
B = 0.5: S2 = 0 ' Начальная точка и сумма квадратов
For J = 1 To JJ ' Перебор точек
QTRAP(A, B, N, C) ' Метод трапеций
Print #2, Using " X = ##.### C(X) = ##.#######"; B;
C
S2 = S2 + (C
-
CT(J))^2 ' Добавка квадрата отклонения
B = B + HB ' Следующая точка
Next J
S2 = S2/N
Print #2, Using " S2 = ##.########"; S2
Print #2,
EndIf
If M = 3 Then
Print #2, " Integration by the Method of Parabol
ic"
B = 0.5: S2 = 0 ' Начальная точка и сумма квадратов
For J = 1 To JJ ' Перебор точек
QPARB(A, B, N, C) ' Метод парабол 13
Print #2, Using " X = ##.### C(X) = ##.#######"; B; C
S2 = S2 + (C
-
CT(J))^2 ' Добавка квадрата отклон
ения
B = B + HB ' Следующая точка
Next J
S2 = S2/N
Print #2, Using " S2 = ##.########"; S2
Print #2,
EndIf
Next M
Close #2
Sleep
' Результаты
вычислений
:
' Integration by the Method of Rectangles
' X = 0.500 C(X) = 0.49244
18
' X = 1.000 C(X) = 0.7819479
' X = 1.500 C(X) = 0.4426246
' X = 2.000 C(X) = 0.4879028
' S2 = 0.00000141
' Integration by the Method of Trapezoids
' X = 0.500 C(X) = 0.4921488
' X = 1.000 C(X) = 0.7757922
' X = 1.500 C(X) = 0.450
5391
' X = 2.000 C(X) = 0.4886461
' S2 = 0.00000561
' Integration by the Method of Parabolic
' X = 0.500 C(X) = 0.4923432
' X = 1.000 C(X) = 0.7799349
' X = 1.500 C(X) = 0.4453544
' X = 2.000 C(X) = 0.4848615
' S2 = 0.00000144
Метод Ромберга
заключается в последовательном уточнении значения интеграла при кратном увеличении числа разбиений. В качестве базовой может быть взята формула трапеций с равномерным шагом H. Если последовательно уменьшать шаг в два раза, получим рекуррентн
ое соотношение для расчета. Теорию метода легко найти в соответствующей литературе.
' P R O G R A M "Quadra2"
' 20.02.2013
'
-------------------------------------------------------------
' Вычислен
ие интеграла методом Ромберга
'
-------------------------------------------------------------
' Ординарная переработка алгоритма "60б" из справочного
' пособия: Библиотека алгоритмов 51б -
100б, М., Соврадио,
' 1976, 136 с., ил.
' В программ
е обозначено:
' X -
переменная интегрирования;
' A -
нижняя граница интегрирования (фиксирована);
' B -
верхняя граница интегрирования (подвижна);
' N -
количество вычисляемых точек функции;
14
' H -
шаг интегрирования;
' S -
вычисленное значение интеграла;
' S2 -
выборочная дисперсия.
#Lang "FB" ' режим совместимости с FreeBASIC
' Вычисление значения подынтегральной функции:
Function Fun(X As Single) As Single
Const As Single Pi = 3.14159265
' X -
аргумент
подынтегральной
функции
Fun = Cos(Pi*X*X/2)
End Function
'
' Подпрограмма интегрирования методом Ромберга
Sub QROMB(A As Single, B As Single, N As Integer, ByRef C As Single)
Dim As Single T(N+1), D, X, H, S
Dim As Integer K, I, J
Dim As LongIn
t M
'
D = B -
A ' Ширина интервала
T(1) = (FUN(A) + FUN(B))/2
K = 1
For I = 1 To N
S = 0: K = 2*K
H = D/K
For J = 1 To K Step 2
X = A + H*J
S = S + FUN(X)
Next J
T(I+1) = (2*S/K + T(I))/2
M = 1
For J = I To 1 Step -
1
M = 4*M ' д
ля типа Integer возможно переполнение!
T(J) = T(J+1) + (T(J+1) -
T(J))/(M -
1)
Next J
Next I
C = T(1)*D
End Sub
'
' Основной модуль программы
Dim As Integer N, J, JJ
Dim As Single A, B, HB, C, CT(4), S2
' Табличные значения функции Френеля:
' C(X) X
Data 0.492344225871 ' 0.5
Data 0.779893400377 ' 1.0
Data 0.445261176040 ' 1.5
Data 0.488253406075 ' 2.0
' Ввод табличных значений
JJ = 4
For J = 1 To JJ
Read CT(J)
Next J
15
N = 8: A = 0: HB = 0.5 ' Количество точек, нижняя граница, шаг
Open "
CON" For Output As #2
'Open "Quadra.res" For Output As #2
Print #2,
Print #2, " Integration by the Method of Romberg"
B = 0.5: S2 = 0 ' Верхняя граница и сумма квадратов
For J = 1 To JJ ' Перебор точек
QROMB(A, B, N, C) ' Метод пря
моугольников
Print #2, Using " X = ##.### C(X) = ##.########"; B; C
S2 = S2 + (C
-
CT(J))^2 ' Добавка квадрата отклонения
B = B + HB ' Следующая точка
Next J
S2 = S2/N
Print #2, Using " S2 = ##.#########"; S2
Print #2,
'
Close #2
Sl
eep
' Результаты вычислений по программе:
' Integration by the Method of Romberg
' X = 0.500 C(X) = 0.49234420
' X = 1.000 C(X) = 0.77989322
' X = 1.500 C(X) = 0.44526118
' X = 2.000 C(X) = 0.48825341
' S2 = 0.000000000
Как это понятно из представленных в программе "Quadra1" методов интегрирования Ньютона
-
Котеса: прямоугольников, трапеций, парабол -
определенные интегралы можно вычислить
,
заранее не задавая число разбиений интервала интегрирования, а указывая необходимую точность вычисле
ния. Ниже приведена программа с замысловатой подпрограммой интегрирования методом Симпсона
с оценкой погрешности.
' P R O G R A M "Quadra3"
' 20.02.2013
'
----------------------------------------------
---------------
' Интегрирование методом Симпсона с оценкой погрешности
'
-------------------------------------------------------------
' на примере вычисления косинус
-
функции Френеля С(х)
' Мудров, стр. 154: 5.6. Вычисление интегралов с заданно
й
' точностью. Программа 5.4B
' Метод Симпсона с оценкой погрешности
#Lang "FB" ' Режим совместимости с QBasic
' Подпрограмма
-
функция вычисления подынтегральной функции
Function Fun(X As Single) As Single
Const As Single Pi = 3.14159265
' X -
аргумент
подынтегральной
функции
Fun = Cos(Pi*X*X/2)
End Function
'
16
' Подпрограмма интегрирования методом Симпсона
Sub QSIMP(A As Single, B As Single, E As Single, ByRef C As Single)
Dim As Single E1, X, H, H1, S, S1, S2, S3
Dim As Integer I
E1 = E*
15: H = (B
-
A)/2
X = A+H: S = 2*Fun(X)
S1 = Fun(A)+Fun(B)+S: S = S1+S
Do
S2 = 0: S3 = S: H1 = H: H = H/2
For X = A+H To B Step H1
S2 = S2+Fun(X)
Next X
S2 = 2*S2: S1 = S1+S2: S = S1+S2
Loop While Abs(1
-
S/(2*S3)) > E1
C = S*H/3
End Sub
'
' Основной модуль п
рограммы
Dim As Integer N, J, JJ
Dim As Single A, B, HB, E, C, CT(4), S2
' Табличные значения функции Френеля:
' C(X) X
Data 0.492344225871 ' 0.5
Data 0.779893400377 ' 1.0
Data 0.445261176040 ' 1.5
Data 0.488253406075 ' 2.0
' Ввод табличн
ых значений
JJ = 4
For J = 1 To JJ
Read CT(J)
Next J
E = 0.00001: A = 0: HB = 0.5 ' Погрешность,нижняя граница, шаг
Open "CON" For Output As #2
'Open "Quadra.res" For Output As #2
Print #2,
Print #2, " Integration by the Method of Simpson"
B = 0.5: S2 =
0 ' Верхняя граница и сумма квадратов
For J = 1 To JJ ' Перебор точек
QSIMP(A, B, E, C) ' Метод прямоугольников
Print #2, Using " X = ##.### C(X) = ##.########"; B; C
S2 = S2 + (C
-
CT(J))^2 ' Добавка квадрата отклонения
B
= B + HB ' Следующая точка
Next J
Print #2, Using " S2 = ###.#########"; S2
Print #2,
'
Close #2
Sleep
'
'
'
17
' Результаты вычислений по программе:
' Integration by the Method of Simpson
' X = 0.500 C(X) = 0.49234319
' X = 1.000 C(X) =
0.77989596
' X = 1.500 C(X) = 0.44526127
' X = 2.000 C(X) = 0.48825315
' S2 = 0.000000000
Любопытный прием демонстрируют методы наивысшей алгебраической точности. В них подынтегральная функция F(X) так же, как и в приведенных выше методах Ньютона
-
Котеса, заменяется полиномами различных степеней. Однако в отличие от методов Ньютона
-
Котеса узлы для построения интерполяционного полинома выбираются из условия обеспечения минимальной погрешности интегрирования. Впервые задача построения квадрату
рных формул подобного типа была решена Гауссом
. Для случая
,
когда пределы интегрирования A и B входят в состав узлов квадратурной формулы, ее называют квадратурой замкнутого типа. Линейным преобразованием всякий конечный отрезок [A, B] может быть преобразо
ван в какой
-
либо стандартный отрезок. Чтобы сделать более простым использование свойств симметрии узлов Xk и коэффициентов Ak, за такой стандартный отрезок мы примем [
—
1, 1]. Ортогональную систему многочленов с постоянным весом на этом отрезке образуют мно
гочлены Лежандра. Квадратуры Гаусса
-
Котеса (Гаусса
-
Кристоффеля) используют например, следующие Ортогональные полиномы:
Полиномы
A
B
W(x)
Лежандра
–
1
1
1
Чебышева
–
1
1
1/Sqr(1
-
x2)
Лагерра
0
+oo
exp(
-
x)
Эрмита
-
oo
+oo
exp(
-
x2)
Оказалось, что вне завис
имости от расположения на отрезке узлов квадратурной формулы она является точной, если подынтегральной функцией является любой многочлен, степень которого на единицу меньше числа этих узлов. Но если есть возможность выбора расположения на этом отрезке неко
торого фиксированного числа N узлов, то ее целесообразно использовать для построения формулы, точной для многочленов возможно более высокой степени. Об этом можно почитать в доступной книге:
Крылов В.И., Шульгина Л.Т. Справочная книга по численному интегри
рованию. M., "Наука", 1966, 372 с., ил.
Ниже представлена программа на языке FreeBASIC, в которой реализован метод Гаусса с шестью узлами Лежандра
.
' P R O G R A M "Quadra4"
' 20.02.2013
'
------------
-------------------------------------------------
' Интегрирование методом Гаусса с шестью узлами
'
-------------------------------------------------------------
' на примере вычисления косинус
-
функции Френеля С(х)
' Мудров, стр. 162: 5.8. Метод
ы наивысшей алгебраической
' точности -
программа 5.7.
18
' Подпрограмма
-
функция вычисления подынтегральной функции
Function F(X As Single) As Single
Const As Single Pi = 3.14159265
' X -
аргумент
подынтегральной
функции
F = Cos(Pi*X*X/2)
End Functi
on
' Подпрограмма метода Гаусса
Sub GAUSS(N As Integer, K As Integer, A As Single, B As Single, X() As Single, C() As Single, ByRef S As Single)
Dim As Single H, A1, B1, R, Q, S1
Dim As Integer I, J
H = (B
-
A)/N
A1 = A
B1 = A + H
S = 0
For I = 1 To N
R = (B1 + A1)/2
Q = (B1 -
A1)/2
S1 = 0
For J = 1 To K
S1 = S1 + C(J)*F(R+Q*X(J))
Next J
S = S+S1
A1 = B1
B1 = B1 + H
Next I
S = S*Q
End Sub
'
' Основная программа
Dim As Integer N
Dim As Single A, B, S
N = 6
Dim As Single X(N), C
(N)
X(1) = 0.932469
X(2) = 0.661209
X(3) = 0.238619
X(4) = -
X(1): X(5) = -
X(2): X(6) = -
X(3)
C(1) = 0.171324
C(2) = 0.360762
C(3) = 0.467914
C(4) = C(1): C(5) = C(2): C(6) = C(3)
A = 0
B = 2
GAUSS(N, 6, A, B, X(), C(), S)
Print Using " S = ###.#########"; S
Sleep
' Результаты вычислений по программе:
' 0.5 0.492344225871
' 1.0 0.779893400377
' 1.5 0.445261176040
' 2.0 0.488253406075
19
3.2.
Вычисление несобственных интегралов.
Несобственные интегралы бывают двух видов:
1
.
Несобственные интегралы с бес
конечным пределом (ами) интегрирования.
Такие несобственные интегралы обычно называют несобственным интегралом первого рода. Для вычисления несобственных интегралов с бесконечными пределами применимы, например, квадратурные формулы Гаусса
-
Кристоффеля, узлы
и веса которых определяются в зависимости от вида весовой функции...
2.
Несобственные интегралы от неограниченных функций.
Такие несобственные интегралы обычно называют несобственными интегралами второго рода. Несобственные интегралы с конечными предела
ми интегрирования, но с подынтегральной функцией, обращающейся в бесконечность в отдельных точках интервала, вычисляют методами аддитивного или мультипликативного выделения особенностей, а также построением нестандартных квадратурных формул.
Известно неск
олько приемов вычисления разных типов несобственных интегралов. Иногда удается заменой переменных перейти от интегралов с бесконечными пределами к интегралам с конечными пределами. Если подынтегральная функция после преобразования останется конечной на нов
ом интервале, то для интегрирования можно использовать методы и программы, рассмотренные выше. В противном случае целесообразно использовать метод, аналогичный методу наивысшей алгебраической точности, но с использованием квадратурных формул Гаусса
-
Кристо
ффеля, но для других ортогональных полиномов.
Для бесконечного интервала интегрирования предлагается квадратура Эрмита с пятью узлами. Шарль Эрмит
(Charles Hermite; 1822 —
1901) —
французский математик, признанный лидер математиков Франции во второй полов
ине XIX века.
' P R O G R A M "Quadra5"
' 20.02.2013
'
-------------------------------------------------------------
' Интегрирование методом Гаусса
-
Эрмита с пятью узлами
'
-------------------------
------------------------------------
' на примере вычисления несобственного интеграла Cos(X)
' Мудров, стр. 168: 5.9. Несобственные интегралы -
' программа 5.8.
' Подпрограмма
-
функция вычисления подынтегральной функции
Function F(X As Single) As Single
Dim As Single A
A = 1
' X -
аргумент подынтегральной функции
F = Cos(X) ' подразумевается *Exp(
-
A*X^2)
End Function
' Подпрограмма метода Гаусса
Sub HERMIT(N As Integer, X() As Single, C() As Single, ByRef S As Single)
Dim As Integer I
20
S = 0
For I = 1 To N
S = S + C(I)*F(X(I))
Next I
End Sub
'
' Основной модуль программы
Dim As Integer N
Dim As Single Pi = 3.14159265
N = 5
Dim As Single X(N), C(N), S, ST, A = 1
X(1) = -
2.020182870456
X(2) = -
0.958572464614
X(3) = 0.000000000000
X(4) = -
X(2)
X(5) = -
X(1)
'
C(1) = 0.019953242059
C(2) = 0.393619323152
C(3) = 0.945308720483
C(4) = C(2)
C(5) = C(1)
'
Open "CON" For Output As #2
'Open "Quadra.res" For Output As #2
Print #2,
Print #2, " Quadrature formula of Gauss"
Print #2, " Orthogonal polyno
mials Ermita"
HERMIT(N, X(), C(), S)
ST = Sqr(Pi)*Exp(
-
0.25) ' истинное значение интеграла
Print #2, Using " S = ###.######### ST = ###.#########"; S; ST
Close #2
Sleep
' Результаты вычислений по программе:
' Quadrature formula of Gauss
' Orthogonal poly
nomials Ermita
' S = 1.380390167 ST = 1.380388498
'
Для полубесконечного интервала интегрирования предлагается квадратура Лагерра с пятью узлами. Эдмон Никола Лагерр
(Edmond Nicolas Laguerre; 1834 —
1886) —
французский математик, член Парижской акаде
мии наук (1885), труды по геометрии, комплексному анализу, исследовал ортогональные многочлены.
' P R O G R A M "Quadra6"
' 20.02.2013
'
-------------------------------------------------------------
' Интегрирование методом Гаусса
-
Лагерра с пятью узлами
'
-------------------------------------------------------------
' на примере вычисления несобственного интеграла Cos(X)
21
' Мудров, стр. 168: 5.9. Несобственные интегралы -
' программа 5.8.
' Подпрограмма
-
функция вычисления подынтегральной функции
Function F(X As Single) As Single
Dim As Single A
A = 1
' X -
аргумент подынтегральной функции
F = Cos(X) ' подразумевается *Exp(
-
A*X)
End Function
' Подпрограмма метода Гаусса
Sub LAGUER(N A
s Integer, X() As Single, C() As Single, ByRef S As Single)
Dim As Integer I
S = 0
For I = 1 To N
S = S + C(I)*F(X(I))
Next I
End Sub
'
' Основной модуль программы
Dim As Integer N
Dim As Single Pi = 3.14159265
N = 5
Dim As Single X(N), C(N), S, ST, A
, B
X(1) = 0.2635603197
X(2) = 1.4134030591
X(3) = 3.5964257710
X(4) = 7.0858100059
X(5) = 12.640800844
C(1) = 0.5217556106
C(2) = 0.3986668111
C(3) = 0.0759424497
C(4) = 0.0036117587
C(5) = 0.0000233700
'
Open "CON" For Output As #2
'Open "Quadra.res" For
Output As #2
A = 1: B = 1 ' параметры функции
Print #2,
Print #2, " Quadrature formula of Gauss"
Print #2, " Orthogonal polynomials Laguerre"
LAGUER(N, X(), C(), S)
ST = A/(A^2+B^2) ' истинное значение интеграла
Print #2, Using " S = ###.######### ST =
###.#########"; S; ST
Close #2
Sleep
' Результаты вычислений по программе:
' Quadrature formula of Gauss
' Orthogonal polynomials Laguerre
' S = 0.500538528 ST = 0.500000000
22
Во время учебы в институте никак не мог понять, почему интеграл от функции COS(X) на бесконечном интервале является величиной конечной
. Как
-
то не слишком отчетливо было сказано о том, что это интеграл от "функции с весом"... Чтобы окончательно убедиться в этом, ниже приведена "итоговая" программа на языке FreeBASIC для вычисления
интегралов понятным методом Симпсона.
' P R O G R A M "Quadra7"
' 20.03.2013
'
-------------------------------------------------------------
' Проверка интегрирования методом Симпсона
'
-----------
--------------------------------------------------
' Используется программа "Quadra3" от 20.02.2013 -
' Метод Симпсона с оценкой погрешности для проверки
' вычислений несобственных интералов первого рода:
' -
Sub HERMIT(N, X(), C(), S) см
. "Quadra5"
' -
Sub LAGUER(N, X(), C(), S) см
. "Quadra6"
#Lang "FB" ' Режим совместимости с QBasic
' Подпрограмма
-
функция вычисления подынтегральной функции
Function Fun(X As Single, WF As Integer) As Single
Dim As Single A, B, WX
' X -
аргумент под
ынтегральной функции;
' WF -
весовая
функция
(Weight Function);
A = 1: B = 1 ' параметры для последующих экспериментов
If WF = 1 Then WX = Exp(
-
A*X^2) ' весовая функция Эрмита
If WF = 2 Then WX = Exp(
-
A*X) ' весовая функция Лагерра
Fun = WX*Cos(X) ' подынтегральная функция
End Function
'
' Подпрограмма интегрирования методом Симпсона
Sub QSIMP(A As Single, B As Single, E As Single, ByRef C As Single, WF As Integer)
Dim As Single E1, X, H, H1, S, S1, S2, S3
Dim As Integer I
E1 = E*15: H = (B
-
A)/2
X = A+H: S = 2*Fun(X,WF)
S1 = Fun(A,WF)+Fun(B,WF)+S: S = S1+S
Do
S2 = 0: S3 = S: H1 = H: H = H/2
For X = A+H To B Step H1
S2 = S2+Fun(X,WF)
Next X
S2 = 2*S2: S1 = S1+S2: S = S1+S2
Loop While Abs(1
-
S/(2*S3)) > E1
C = S*H/3
End Sub
'
' Основной мо
дуль программы
Dim As Integer N, WF, WFM
WFM = 2 ' два типа полиномов
Dim As Single Xmin(WFM), Xmax(WFM), E, C
E = 0.0000001 ' Погрешность
вычисления
23
Open "CON" For Output As #2
'Open "Quadra.res" For Output As #2
Print #2,
Print #2, " Integration by the Method of Simpson"
Xmin(1) = -
5: Xmax(1) = 5 ' Верхняя и нижняя границы
Xmin(2) = 0: Xmax(2) = 15 ' Верхняя и нижняя границы
For WF = 1 To WFM
QSIMP(Xmin(WF), Xmax(WF), E, C, WF) ' Обращение к подпрограмме
Print #2, Using " WF
= ## Xmin = ##.### Xmax = ##.### C(WF) = ##.########"; WF; Xmin(WF); Xmax(WF); C
Next WF
Print #2,
'
Close #2
Sleep
' Результаты вычислений по программе:
' Integration by the Method of Simpson
' WF = 1 Xmin = -
5.000 Xmax = 5.000 C(WF) = 1.3803885
0
' WF = 2 Xmin = 0.000 Xmax = 15.000 C(WF) = 0.50000024
Будем живы, нужно обязательно рассмотреть тему: "
Классические ортогональные полиномы"
–
она сулит множество интересного...
Конечно, учитывая возросшее быстродействие компьютеров, вычисление простых и кратных интегралов можно осуществить повторным применением алгоритмов и программ, рассмотренных в предыдущих разделах, но весьма соблазнительно использовать всевозможные варианты метода статистических испытаний (метода Монте
-
Карло), а с этим мето
дом мы уже знакомы: "Фрагмент 2. Все еще впереди" содержал программу №4 для вычисления числа Pi методом Монте
-
Карло. С различными вариантами методов Монте
-
Карло, их теоретическим обоснованием можно ознакомиться по книге:
Ермаков С.М., Михайлов Г.А. Статист
ическое моделирование. Том 1: Наука, 1982. 296 с.
Можно, при большом желании
,
познакомиться с программами, приведенными Научно
-
образовательным Интернет
-
ресурсом НИВЦ МГУ по численному анализу:
http://num
-
anal.srcc.msu.ru/lib_na/cat/cat7.htm
Раздел "Числен
ное интегрирование".
Однократные интегралы (без контроля точности).
Однократные интегралы (с контролем точности).
Многократные интегралы.
Интегралы специального вида.
Интегралы с особенностями и несобственные интегралы.
Мне, например, понравилась подпрогр
амма QSL6R: "Вычисление определенного интеграла по формулам Лобатто 11
-
ой степени
точности для больших отрезков интегрирования от функций с локализованной особенностью":
http://www.srcc.msu.su/num_anal/lib_na/cat/q/qsl6r.htm
но уж больно она мудреная...
24
На этом позвольте закруглиться.
До новых встреч!
Пишите: eugene_r@mail.ru
Документ
Категория
Информатика
Просмотров
189
Размер файла
575 Кб
Теги
интеграл, алгоритм, программа, freebasic
1/--страниц
Пожаловаться на содержимое документа