Fedorenko-RP-Vvedenie-v-vychislitelnuyu-fiziku (810773), страница 82
Текст из файла (страница 82)
Тогда Л(~ У) = иао/с + тгьв А~нс + иа ~д +|+ м|,| А|н, | (1о) где ив Ф= Ь з(1Ь+ Ь вЂ” хНтЬ+ Ь вЂ” у) и т.д. Элеменг «матрицы» 1 (1, 1; х, у), очевидно, вычисляется так: по точке (х, у) нужно найти индексы 1, т и коэффициенты интерполяции, после чего определюотся элементы матрицы 1'.
1'(1+ а, т + Р; х, у) = и, 8, а, Р = О, 1. Все остальные элементы 1 равны нулю. Итак, операция раздачи заряда аы находящегося в точке (Х, г ) состоит в следующем. По значениям Хы Уь находится ячейка сетки (1, т), в которой эта точка находится, вычисляются соответствующие коэффициенты интерполяции и 8, и доли заряда а«|г раздаются в узлы (1+ а, т + Р): Р|+а,а»а' Р|»а,е»В+ |4 ар Если бы использовалась более точная интерполяция, мы столкнулись бы с таким «нефизичным» явлением, кзк раздача отрицательной доли заряда в какой-то узел. Видимо, поэтому на практике ограничиваются простейшим кусочно-билинейным оператором интерполяции, когда все и 8 |в О и ~ч' |г 8— - 1.
Таким образом конструируются схемы интегрирования, обеспечивающие непрерывность плотности заряда р, по положениям частиц (Хы г'„) и сохранение импульса. Быстрое дискретное преобразование Фурье. Имеется некоторая сеточная функция 1(Ь) (Ь = О, 1, ..., Ф вЂ” 1). Она может быть представлена разложением в сумму Фурье: У-| 1(Ь)= ' 7 1(л) ., = '", Ь=О,1,...,Ы 1, (20) и 0 з 241 ПРИБЛИЖЕННОЕ ИНТЕГРИРОВАНИЕ РРАВНЕНИЯ ЕЛАСОБА ззо Коэффициенты Фурье А(п) вычисляются по аналогичным формулам; И-1 А(п) =д~~ У(/с) 1Р "", и =О, 1, ..., Дà — 1.
(21) «-о Формулы (21), (20) представляют собой прямое и обратное преобразования. Выполнение каждого из них требует О(ЛГЕ) операций, если программировать непосредственно формулы (20), (21). Заметим, что нужно аккуратно вычислять степени 1Р, чтобы избежать лишних операций. Ниже это будет сделано. Однако наиболее существенное сокращение обьема вычислений связано с тем, что при вычислении сумм многие слагаемые и даже группы слагаемых вычисляются неоднократно, и этого повторения можно избежать. На этом основан алгоритм РРТ (Разс Роппе Тгапз1огщ). Рассмотрим реализацию (20), опуская, для простоты, множитель 1/т'М.
Заметное сокращение операций происходит тогда, когда чис- лО ЛГ можно разложить иа множители. Пусть лс' = МттГ1. Представим lс и и в виде 'с ХФ1+й1 и 9+ п1М (гг) Здесь И=0,1,...,М вЂ” 1;е1=0,1,...,ЛГ1 — 1; 9=0,1,..., М вЂ” 1; и, = О, 1, ..., М1 — 1. Показатель степени еп можно записать так: ел = й(9+ п,М) = ет+ 1сп,М= ЕР+ псМ(хлс1 + Гсс) = = ЕР+ п1ИМХ1 + ксп1М. Поскольку ими = 1ои = 1, слагаемое псхМПГ1 можно из хп убрать и формула (20) примет вид И-1 и-! Д/с) = ~ ~А(п) ме сое~" и = ~ А(п) 1Р~" РРЕ~~"ч и' = тгм. (20а) -о Введем двухиндексную нумерацию, используя одно и то же имя для функций с разным числом целочисленных аргументов: У(х, 11) Бяу(х1У1+ йс), А(Р, и,) юА(9+ п,М). (23) Сумму по и ° О, 1, ..., Лà — ! представим как двойную: И-1 И,-1 (20б) У(х, й1) = ~~; ~ А(Р, п1) н"" и~ "и я о я о с где и =О, 1, ..., М вЂ” 1; Е1 =О, 1, ..., Ф1 — 1. Внутренняя сумма в (20б) вычисляется при постоянных Р и е, поэтому множитель НРР' пгмвлнжвнныв методы вычислительной ьнзнкн 1ч.
и можно вынести за скобки: и-~ ! А(т, л~) мф~1м1 м ~". м 0 1 и-! У(х, м,) = ~ «=0 (20в) Обозначим внутренние суммы, зависящие лишь от т и й~.' и-~ 1 А,(т, Ус,) ='~ А(т, л,) мьР, (24) м О 1 т=0,1,...,М вЂ” 1, /с,=0,!,...,М,— 1. Тогда (25) У(х, А,) = ~~, А,(т, А,) нь', х=0,1,...,М вЂ” 1, 1,=0,1,...,М,— 1, ~=хи,+й,. Итак, преобразование (20) распалось на два последовательных. Сначала исходные величины А преобразуются в А, по формуле (24); это стоит СММьМ, операций (постоянную С вычислим ниже: С = О(1)). Затем за СМММ, операций выполняется преобразование (25). Итого мы имеем С1т(М+ Ф,) операций вместо СИМА,. При больших М и МжФ, ~(Ф это уже заметный выигрыш.
Если число Ж, может быть разложено на множители, этот прием может быть использован снова. Будем понимать (24) как М преобразований типа (20). В самом деле, единственное свойство основания н в (20), которое было использовано, есть соотношение мн=1.Но и м~» =м"'и =1. Пусть М=М, — простое число, а М, = Мхах Тогда при каждом т = О, 1, ..., М1 — 1 нужно вычислить Ф, коэффициентов А,(т, /с,), а зто можно сделать за СЖ,(М + Ф ) операций. Таким образом, мы имеем СФ,М,(Мз+Мз) операций. Вместе с СМ,М,М, операциями на вычисление (25) получаем СФ(М, + М +Ж ) операций.
И т.д. Наибольший эффект достигается в случае, когда Ф разлагается на множители малой величины (2, 3, например). Наиболее популярен алгоритм ггТ при Ф= 2г. В этом случае число операций сводится к СФ(2р) = 2СФ 1ок К. Для оценки постоянной С приведем текст программы, реализующей, напри- о 241 пгивлижаинох интеггиговлиие зтланквяя вллсовл 39! мер, преобразования (25). Используем нестрогую версию языка РОСТКА!л(: "О = "' ' » "1 = ! л бе!к=О,М вЂ” ! й до 2 к1= О, А11 — 1 3 11з= ! 4!о 3 т = О, М вЂ” 1 $ Дх, ~с1) = У(х, й1) + Аз(т, й1)»оз "З = "З""г*" 1 Ог = Ог'1» "1 = "1»ОО Заметим, что и""" = и4"'9+4 >", а величины ви иг, оз в циклах по х, к„т принимают значения „, Лз, „,„4, („»„)»,„14,+»Л1Р— о — г — з — 1 г Решение уравнения Пуассона методом Фурье. Разиостиое уравнение Пуассона (7) легко решается методом Фурье.
Обозначим = 4хр,"+ 'йг. Сначала делаем преобразование по первому индексу. Получаем прсдставление У1 в виде у =~Гл л Это А1 быстрых преобразований (при каждом лз отдельно), которые выполняются за СА1» !ойг М операций (при М = 2»). Затем при каждом Л выполняем преобразование по индексу р. Получаем выражение У = а С' ~' Рл» 1»1л+ г, а=!уз/У, л юи о Зто спзит еще СМг 1ойг М операций. В з 14 было показано, что сеточные функции згзл+» являются собственными функциями разностного оператора Лапласа. Соответствующие собственные значения суть Л = — — 1 з)п — + зйп 4( .
гЛ» . гилЛ Ог ! 2Л1 241) Из (7) получаем коэффициенты Фурье искомой функшзи ~Р: Фл» = Р"" о/Лл . После этого дважды делаем обратное преобразование и находим значения в узлах сетки: р,„,=а~~ ~Х,Фл л» Итак, решение (7) стоит 4СА12!охг М операций. пеивлижвнныа методы вычислительной физики 392 1Ч. и й 2з. Некорректные задачи н нх прнбпнженное решение Почти все задачи, с которыми мы имеем дело, можно записать в общей форме: 22(и) =/, (1) где и — искомая функция из некоторого функционального пространства У, / — известная «правая часть» уравнения, принадлежащая пространству Р, а Я вЂ” оператор (вообще говоря, нелинейный) из У в Р.
Обычно правая часть известна нам не точно, а с некоторой погрешностью Ь, т.е. «реальной действительности» соответствует «истинная» функция /, а решается задача с правой частью /, причем й/ — /й„-~ Ь (где й 'й — норма в пространстве Р), Естественно, возникает вопрос: в какой мере погрешность в / сказывается на решении, т,е. каково отличие и от й, где й — решение «точной» задачи Я(й) =/? Определение 1. Задача (1) называется поставленной корректно, если из й/ — /йг м Ь следует зй — ийи ц з, где «может быть сделано сколь угодно малым за счет достаточно малого Ь. Обычно считается, что физические задачи приводят к корректно поставленным математическим и только последние подлежат приближенному решению. Наоборот, задача считается поставленной некорректно, если при сколь угодно малом Ь и сколь угодно большом Л найдется функция /, такая, что при ~~/ — /й' ж Ь оказывается йй — ий > Л.
Есть ли смысл находить решение некорректной задачи, если сколь угодно малые погрешности в правой части (а они всегда есть) приводят к большим погрешностям в решении? Заметим, что характеристика задачи как корректной или некорректной зависит от выбора норм й ~~и и й й ., и задача, некорректная при одном выборе норм, может оказаться корректной при другом. Однако зтот формальный способ «исправления» задачи обычно не проходит, так как выбор норм не является совершенно произвольным и должен правильно отражать суть дела в содержательной постановке задачи.
В дальнейшем мы поясним эти вещи на конкретном примере. Рассмотрим классический пример (Адамара) некорректной задачи — так называемую обратную задачу теплопроводности. Ищется функция и(г, х), удовлетворяющая уравнению (2) з 251 НЕКОРРЕКТНЫЕ ЗАДАЧИ И ИХ ПРИБЛИЖЕННОЕ РЕШЕНИЕ зоз в прямоугольной области 0 И х ~ 1, 0 а:.
Г М Т, с начальными данными и(0, х) = /(х) и краевыми условиями и(0 0) = и(0 1) = О. Пусть й(с, х) — решение этой задачи. Возмутим «правую часть» (входную информацию), т.е. рассмотрим ту же задачу, но с начальными данными и(0, х) = /(х) + Ь яп (асскх), где Ь вЂ” очень малое число. Выпишем явно решение новой задачи: и(а, х) = и(0 х) + Ь ег ' ' яп (екх). При достаточно большой частоте к возмущение решения может стать сколь угодно большим при сколь угодно малом возмущении начальных данных. Покажем, как сделать эту задачу корректной, выбирая нормы подходюцим образом.