О законе сохранения замолвите слово
Банальные, в общем-то, рассуждения о свойствах численных решений дифуров.
Пришлось на днях решать до боли простое и известное трехмерное сферически-симметричное уравнение диффузии с источником
&partt f =D&partR(R2&partR(f))/R2 +Q(t)&delta(R)
Гран. условия есть конечно, но тоже самые простые, f(R=&infin, t )=0 , f(R,0)=0 , &partRf(R=0,t)=0 .
Казалось-бы, что сложного -- от схемы многого не хочется, пусть будет самая простая, явная.
Заменяем производные правыми разностными. Вместе с условием устойчивости решения (что-то вроде dt < dr2/D) получаем схему первого порядка по пространству и времени.
Пишем схему, получаем решение, радуемся. Но преждевременно. Легко посчитать количество инжектированных за время T частиц: &int Q(t)dt, а с другой стороны это равно 4&pi &int f(R,T)R2dR.
Давайте теперь посчитаем эти два интеграла-суммы в нашей схеме. Окажется (surprise,surprise!), что они и близко не равны.
Hint: считать удобнее всего разницу &Sigmar ( f(t+dt,r)-f(t,r) )wr , где wr -- некоторые веса, wr ~ r2 -- разницу количества частиц на шаге t и t+dt
Еще один сюрприз -- для одномерного уравнения для аналогичной схемы, эти суммы совпадают. Может из-за того, что в одномерном случае схема оказывается более точной по пространству (2-го порядка)? Ан нет, даже если заменить в 3-х мерном случае схему на более точную (читай правую производную -- на центральную) -- не помогает.
Интересно?
Решение оказывается простым до безобразия -- хотим сохранения количества частиц, надо делать его самим.
Переписываем дифур, интегрируя левую и правую части по R,t по клеточке (R-dr,t ; R,t+dt) (левый-нижний, правый-верхний углы). Заменяем интегралы суммами, получаем разностную схему похожую на ту, что была с самого начала, плюс поправка порядка (f(t,R)-f(t,R-dr))/R2 (при том, что &partrf~1/dr).
Как же получается так, что начальная схема не сохраняет количество частиц? Оказывается просто. Схема при уменьшении шага сходится к правильному решению, но идет, например, везде выше его. Тогда в каждой точке разница между решениями маленькая, а между интегралами решений (читай -- закон сохранения) может быть большой.
Схемы, которые сохраняют интегралы уравнения люди называют консервативными. Простая замена производных на сколь-угодно точные разностные, похоже не гарантирует консервативности схемы :-(
Вот почему нам обо всем этом не рассказывали на каких-нибудь численных методах, а? Или на пары надо было больше ходить?
Да, чуть не забыл, я же уже билеты домой взял. Так что перемещения:
Дублин -->Киев : 22.06
Киев -->Дублин : 15.08
В промежутке тоже надеюсь куда-нибудь выбраться, хоть ненадолго :-)
Пришлось на днях решать до боли простое и известное трехмерное сферически-симметричное уравнение диффузии с источником
&partt f =D&partR(R2&partR(f))/R2 +Q(t)&delta(R)
Гран. условия есть конечно, но тоже самые простые, f(R=&infin, t )=0 , f(R,0)=0 , &partRf(R=0,t)=0 .
Казалось-бы, что сложного -- от схемы многого не хочется, пусть будет самая простая, явная.
Заменяем производные правыми разностными. Вместе с условием устойчивости решения (что-то вроде dt < dr2/D) получаем схему первого порядка по пространству и времени.
Пишем схему, получаем решение, радуемся. Но преждевременно. Легко посчитать количество инжектированных за время T частиц: &int Q(t)dt, а с другой стороны это равно 4&pi &int f(R,T)R2dR.
Давайте теперь посчитаем эти два интеграла-суммы в нашей схеме. Окажется (surprise,surprise!), что они и близко не равны.
Hint: считать удобнее всего разницу &Sigmar ( f(t+dt,r)-f(t,r) )wr , где wr -- некоторые веса, wr ~ r2 -- разницу количества частиц на шаге t и t+dt
Еще один сюрприз -- для одномерного уравнения для аналогичной схемы, эти суммы совпадают. Может из-за того, что в одномерном случае схема оказывается более точной по пространству (2-го порядка)? Ан нет, даже если заменить в 3-х мерном случае схему на более точную (читай правую производную -- на центральную) -- не помогает.
Интересно?
Решение оказывается простым до безобразия -- хотим сохранения количества частиц, надо делать его самим.
Переписываем дифур, интегрируя левую и правую части по R,t по клеточке (R-dr,t ; R,t+dt) (левый-нижний, правый-верхний углы). Заменяем интегралы суммами, получаем разностную схему похожую на ту, что была с самого начала, плюс поправка порядка (f(t,R)-f(t,R-dr))/R2 (при том, что &partrf~1/dr).
Как же получается так, что начальная схема не сохраняет количество частиц? Оказывается просто. Схема при уменьшении шага сходится к правильному решению, но идет, например, везде выше его. Тогда в каждой точке разница между решениями маленькая, а между интегралами решений (читай -- закон сохранения) может быть большой.
Схемы, которые сохраняют интегралы уравнения люди называют консервативными. Простая замена производных на сколь-угодно точные разностные, похоже не гарантирует консервативности схемы :-(
Вот почему нам обо всем этом не рассказывали на каких-нибудь численных методах, а? Или на пары надо было больше ходить?
Да, чуть не забыл, я же уже билеты домой взял. Так что перемещения:
Дублин -->Киев : 22.06
Киев -->Дублин : 15.08
В промежутке тоже надеюсь куда-нибудь выбраться, хоть ненадолго :-)