LINUX.ORG.RU
ФорумTalks

Олимпиада LOR

 олимпиада lor, ,


0

6

Друзья!

Почитав тему про определение эффективности ВУЗов, да и просто ради развлечения, предлагаю устроить небольшое состязание по численным методам.

Приглашаются все желающие, способные решить задачу, требующую применения численных методов. Уровень методов  — от простейших, до весьма навороченных — позволяет по разному продвинуться в точности решения. Правила просты: решите задачу, написав программу на любимом языке, и запостите сюда три числа в качестве ответа с тем числом значащих цифр, которые считаете верными (сошедшимися). Победителем будет считаться тот, кто первым запостит наиболее точный ответ. Предельной точностью считается точность double — 14-15 значащих цифр. Программу не проверяю — только финальный ответ, но победитель может (если захочет) запостить ее сюда — это одобряется и приветствуется. Сама задача идет ниже.

Задача. Решение уравнение Пуассона-Больцмана для функции распределения противоионов в растворе вокруг центральной (коллоидной) частицы.

В центре сферической полости радиуса R=5σ находится коллоидная частица радиуса σ с зарядом Z = 20. Ее противоионы имеют пренебрежимо малый размер и могут быть где угодно в растворе от σ < r < R. Найти функцию распределения противоионов g(r) в растворе, предполагая, что они подчиняются уравнению Пуассона-Больцмана. А именно:

1. Функция распределения g(r) пропорциональная больцмановской экспоненте от электростатического потенциала, ощущаемого частицей на расстоянии r:
     g(r) ~ exp(-ξ Ueff(r)),          (1)
где ξ -- длина Бьеррума (параметр задачи, который будет задан ниже). Коэффициент пропорциональности определяется из условия электронейтральности системы ("нормировки")
     \int_σ^R 4πr²g(r) dr = Z.        (2)

2. Электростатический потенциал, в свою очередь, определяется как:
     Ueff(r) = -Qin(r)/r + \int_r^R 4πr g(r) dr,  (3)

где 
     Qin(r) = Z - \int_σ^r 4πr² g(r) dr.          (4)

Решить задачу для случая ξ = 1/2 σ.

Вопросы:

а) Найти функцию распределения g(r). Поскольку сетки у всех могут быть разными, в качестве результата запостить g(R) [просто] и g(σ) [сложно].

б) Найти полную электростатическую энергию системы для полученного равновесного распределения:

E = -\int_σ^R Qin(r) 4πr g(r) dr.

★★★★★
Ответ на: комментарий от unanimous

Ну так тема-то про эффективность ВУЗов, никто школьные задачки тут давать не будет.

Да я не про то, что задачи дожны быть школьные. Я про очередное доказательство того факта, что когда ученый долго варится в своей среде, ему начинает казаться что любой думающий человек должен разбираться в его области. «Ты физик?» — «Ну да, вроде физик» — «И ты не разбираешься в этой... как там Alv ее назвал? физической кинетике?» — «Даже слов таких не слышал» — «Да какой же ты физик после этого?»

А в конце делаются выводы что ЛОР полон школьников, и что здесь только фары сверхсветовых ракет умеют обсуждать. Ну или на худой конец черные дыры.

morse ★★★★★
()
Ответ на: комментарий от AIv
#######################################################################
## This program finds Tc for ferroelectric/paraelectric transition.
## Ferroelectric material has |<x>| > 0 (spontaneous polarization), so
## Tc is determined as the lowest temperature when <x> becomes 0.
#######################################################################

clear;

T      = 0.3865;		# Termodynamic temperature
a = 1; b =0.1; alpha = 1.1;	# U(x) paratemers
Npts   = 128;			# Grid size
Niter  = 1000000;		# Maximum number of iterations
errtol = 1e-10;

## We are working on Clenshaw-Curtis grid. \pm 10 replaces \pm infinity
[x, w] = clencurt(Npts, -10, 10);

f = U = zeros(1, Npts);

## <x> initial guess
xav = 1.0;

tic
for iter=1:Niter
  xav_old = xav;

  U = a*x.^2/2 + b*x.^4/4 - alpha*xav*x; # Potential
  f = exp(-U/T)/(sum(exp(-U/T).*w));	 # Distribution function
  xav = sum(x.*f.*w);			 # <x>

  ## Convergence test
  if ((abs(xav-xav_old)/xav < errtol) || (xav < errtol))
    fprintf("xav converged@iter = %d, xav = %.14g\n", iter, xav);
    break;
  endif

endfor
toc

## Answer will be: 0.3863 < Tc <= 0.3864
unanimous ★★★★★
() автор топика
Ответ на: комментарий от morse

Я про очередное доказательство того факта, что когда ученый долго варится в своей среде, ему начинает казаться что любой думающий человек должен разбираться в его области. «Ты физик?» — «Ну да, вроде физик» — «И ты не разбираешься в этой... как там Alv ее назвал? физической кинетике?» — «Даже слов таких не слышал» — «Да какой же ты физик после этого?»

Я в жизни не имел дело с пьезоэлектриками — как видишь, это не помешало мне решить задачу. В этом-то вся и соль: ни в моей задаче, ни в задаче Alv не требуется глубокого понимания предметной области, скорее общая физическая культура. И моё ощущение, что именно эта культура важна. Есть такая замечательная поговорка: «Образование есть то, что остается, когда всё выученное уже забыто», это именно про ту самую культуру предметной области.

К тому же, участие было свободное — никто и не ожидал от программиста там или электронщика решения, хотя бы путь физики решили бы.

unanimous ★★★★★
() автор топика
Ответ на: комментарий от unanimous

Это какой ЯП? Так то логику я вроде понял.

Аналитическое решение дает Тс=1.1/3 = 0.3(6).

Но в окрестностях крит. температуры все сходится очень долго по понятным причинам, там и по физике время релаксации стремится к бесконечности. У меня аналогичный код давал Тс~4 (чуть меньше) - точка перехода размыта и не очень понятно что считать за Tc. Если max|d<x>/dT| то примерно 3.7-3.8

AIv ★★★★★
()
Ответ на: комментарий от AIv

Это какой ЯП? Так то логику я вроде понял.

Octave, после небольшой обработки grep-ом должен завестись на Matlab'е.

unanimous ★★★★★
() автор топика
Ответ на: комментарий от AIv

clencurt.m генерирует узлы и веса квадратуры Clenshaw-Curtis'а — универсальная квадратура, обеспечивающая экспоненциальную сходимость интегральных сумм для бесконечно-гладких функций и полиномов.

unanimous ★★★★★
() автор топика
Ответ на: комментарий от morse

Человек, не знающий трех законов Ньютона, начал термодинамики, распределений Максвелла и Больцмана и пр. «основ мироздания» (без вывода, просто о чем речь идет хотя бы) - не физик ни разу. Называть он себя при этом может как угодно, хоть папой римским, и диплом у него может быть любого вуза любого цвета.

Ну а когда уравнения записаны - тут уже не важно что это, физика, генетика или лингвистика. Математика тем и хороша, что это универсальный язык науки.

У ТС кстати задача сложнее моей первой, но проще моей второй. Ну если о точности не задумываться... ;-)

AIv ★★★★★
()
Ответ на: комментарий от unanimous

Octave, после небольшой обработки grep-ом должен завестись на Matlab'е.

ни того ни другого нету... у нас плюсы да питон. Ок, я чуть разгребу и на выходных постараюсь Вашу задачу таки решить. А то что ж, зря все было что ли;-)

AIv ★★★★★
()
Ответ на: комментарий от unanimous

В окрестностях Tc область, где все очень медленно сходиться. В итоге для любой заданной точности получается «хвост» в сторону больших Tc температур.

AIv ★★★★★
()
Ответ на: комментарий от AIv

Человек, не знающий трех законов Ньютона, начал термодинамики, распределений Максвелла и Больцмана и пр. «основ мироздания» не физик ни разу.

Во-первых, вот я например уже сто лет как забыл всевозможные распределения, а из начал термодинамики помню только второе, и то потому что оно скорее философское, нежели физическое. Во-вторых, то что ты сейчас написал является прямым подтверждением моих слов: «Как? Ты не разбираешься в моей области?» и далее по тексту. Почему именно статфизика должна являться «основой мироздания»? Мне вот например куда ближе какие-нибудь диаграммы Фейнмана. А ведь сколько еще есть областей, которые с данной темой никак не пересекаются.

morse ★★★★★
()
Ответ на: комментарий от AIv

поломаю голову на досуге. Не удивляйся если ответ на проверку пришлю через пол года :)

true_admin ★★★★★
()
Ответ на: комментарий от morse

В том то и дело, что перечисленное мной не является «моей» областью, а является неким необходимым минимумом знаний о природе. В отличии скажем от диаграмм Фейнмана - я представляю о чем речь, и не более того, таки 99.(9)% тоже «не более того». Никто ж не требует сходу писать цепочку ББГКИ какую нить... и даже у-е Больцмана или даже Власова. Никто не требует рассказывать про приближение Хартри-Фока или какие нить матрицы Паули.

Но уж ЗСЭ, ЗСИ, ЗСМИ, какие то общие принципы СТО и ОТО, что есть потенциал и что есть поле (электрическое и магнитное), как выглядит плоская ЭМ-волна (бог бы с ними с у-я Максвелла), и пр. «базовые» вещи любой физик знать обязан. И неважно относится это к его области или нет, работает он физиком или дворником. Это остается навсегда, это действительно элемент культуры человека, нормально учившегося в каком то околофизическом вузе.

AIv ★★★★★
()
Ответ на: комментарий от AIv

Золотые слова. А теперь бери эти «базовые» вещи — законы сохранения, общие принципы и т.д. и попытайся решить хоть одну из вышеперечисленных задач только ими. Или вообще хоть какую-нибудь задачу требующую точного решения. Быстро окажется, что базовые вещи хороши только если надо на пальцах объяснить как GPS работает (и то не у всех выйдет).

morse ★★★★★
()
Ответ на: комментарий от morse

Для решения «этих» задач вообще ничего кроме знания тупого метода посл итераций и умения интегралы численно считать (от известной аналитической ф-ии) не требуется. Базовые вещи нужны, что не впадать как Eddy_Em в ступор при виде слова «электрик»;-)

Кроме того, понимание физики дела существенно облегчает численное решение и позволяет проверять промежуточные результаты на предмет «здравого смысла». «Чистый» математик любую задачу решит правильнее, но провозится куда дольше, поскольку физик будет срезать углы основываясь на каких то общих соображениях.

AIv ★★★★★
()

Ну как то так... g(sigma)=2.833, g(R)=0.01925, E=-130.7 Считал что sigma=1, размерность исходных у-й не проверял. За 100 итераций ошибка g падает до 10^-15, но влияние размера сетки тоже не проверял. Я был не прав, задача не такая скучная. В частности тупо последовательными итерациями не сходиться, возникает болтанка, почему пока не понял.

#!/usr/bin/python
from math import *

Z, R, xi = 20, 5., .5
N = 1000
#-----------------------------
gL, dr = [ 3*Z/(4*pi)/(R**3-1) ]*N, (R-1.)/(N-1)

for i in range(100):
    Qin, Uout, gL2, E = Z, 4*pi*( sum([ (1+i*dr)*g for i, g in enumerate(gL) ])-.5*(gL[0]+R*gL[-1]) )*dr, [], 0
    for i, g1 in enumerate(gL[:-1]):
        r1 = 1+i*dr; r2 = r1+dr; g2 = gL[i+1] 
        gL2.append(exp( -xi*( -Qin/r1 +Uout ) ))
        Q1 = Qin; Qin -= 2*pi*( r1*r1*g1 + r2*r2*g2 )*dr
        Uout -= 2*pi*( r1*g1 + r2*g2 )*dr
        E -= 2*pi*( Q1*r1*g1 + Qin*r2*g2 )*dr
    gL2.append(exp( xi*Qin/R ))
    _norm = Z/( 4*pi*( sum([ (1+i*dr)**2*g for i, g in enumerate(gL2) ])-.5*(gL2[0]+R*R*gL2[-1]) )*dr )
#    gL = [ g*_norm for g in gL2 ]
    gL = [ .5*( g1+g2*_norm ) for g1, g2 in zip(gL,gL2) ]
    acc1 = max([ abs(g1-g2*_norm) for g1, g2 in zip(gL,gL2) ])
    print gL[0], gL[-1], _norm, E, acc1, Uout, Qin
#for i, g in enumerate(gL): print i, g
AIv ★★★★★
()
Ответ на: комментарий от AIv

Да, верно. И осцилляции там есть. Но на равномерной сетке трудно получить хорошую точность, хотя вот тебе довольно хорошее значение g(\sigma) удалось сосчитать.

unanimous ★★★★★
() автор топика
Ответ на: комментарий от unanimous

На самом деле, что мне очень нравится в этой задаче — то, что она многоуровневая. На первом уровне нужно ее хотя бы понять. На втором — транслировать в код, пободаться с осцилляциями, получить решение с небольшой точностью. На третьем, последнем, если хорошо подумать и использовать продвинутые подходы можно получить решение с очень высокой точностью. Очень, очень приятная задачка.

Плохо только то, что лишь один человек со всего ЛОРа с нею справился. Но таки за труды спасибо.

unanimous ★★★★★
() автор топика
Ответ на: комментарий от unanimous

Да ладно. Зачем в этой задаче высокая точность то? Три знака ИМНО за глаза... модель то простенькая. Да и g ведет себя нормально, без каких то выбросов и рывков, совсем не очевидно что нужна неравномерная сетка. Меня всегда удивляло, как много трудов народ тратит на получение точного решения в приближенной постановке;-)

ИМНО тут интересней посмотреть на динамику установления решения (хотя наверное вы такими вещами пренебрегаете, считаете что равновесие достигается сразу?) или скажем на вклад во взаимодействие двух частиц со стороны противоионов. Вообще на плазму похоже...

AIv ★★★★★
()
Ответ на: комментарий от AIv

Да ладно. Зачем в этой задаче высокая точность то?

Для проверки того, что из себя представляет решающий :) Умение найти хороший метод решения — часть теста.

unanimous ★★★★★
() автор топика
Ответ на: комментарий от unanimous

Для проверки того, что из себя представляет решающий :)

Отец преподавал в инженерном вузе, ЕМНИП за точность больше трех знаков снижал оценку. Тест хороший, кто выдал весь дабл тот дурак?;-)

AIv ★★★★★
()
Ответ на: комментарий от AIv

ЕМНИП за точность больше трех знаков снижал оценку. Тест хороший, кто выдал весь дабл тот дурак?;-)

Я сказал рапортовать *сошедшиеся* знаки. Их в наивной формулировке получать очень трудно: сходимость линейная (или близкая к таковой), поэтому сетки требуются гигантские. Или правильный метод решения :)

unanimous ★★★★★
() автор топика
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.