Лекция Компютърно моделиране във физиката роля

Етапи в компютърното моделиране

Написването на една компютърна програма до голяма степен повтаря етапите на един експеримент. Не случайно наричат компютърното моделиране още числен експеримент. Написването на програмата и изчистването и от синтактични и логични грешки е аналогично на създаването на експериментална установка. Наличието на полезен сигнал от установката обаче не означава, че експеримента дава коректни резултати. По същият начин фактът, че една програма няма синтактични грешки (т.е. програмата дава някакъв резултат) не означава, че тя работи правилно. Затова анализът на възможните грешки и факторите, които стоят в тяхната основа е неизменен и до голяма степен най-важен етап във всеки физически експеримент. Експерименталните грешки биват случайни и систематични. Първите трябва да бъдат сведени близо до възможния теоретичен минимум, а вторите трябва да бъдат отчетени и доколкото е възможно да бъдат избегнати. В този случай аналогията с компютърното моделиране е пълна: един грешен алгоритъм може да даде недопустими систематични грешки, а неправилно подбрани параметри на програмата (например стъпката във времето или недостатъчна точност) могат да доведат до флуктуации много по-големи от теоретично предсказаните.

Също както при сложен експеримент отделните части на една програма трябва щателно да се тестват преди да се сглобят в едно цяло. След сглобяването на една установка и проверката на получените резултати следва тяхния анализ. Нито един експеримент не се поставя сам за себе си, извън рамките на определени хипотези и резултатите винаги се интерпретират в рамките на даден модел, било за да го отхвърлят или за да го потвърдят. По същия начин се анализират резултатите от числения експеримент. Съпоставянето с теоретичните резултати, получени от опростени модели показва първо доколко валидни са направените приближения и второ, по-важно, позволява да се оценят новите резултати и да се навлезе по-дълбоко в същноста на изучаваните явления.

Въведение в метода на компютърния експеримент (КЕ)

В основата на КЕ лежи определен модел на физическата система, чийто характеристики искаме да определим. Тези характеристики (или наблюдаеми величини) се получават след усредняване по пространството на състоянията на системата. Нека означим хамилтониана на модела на разглежданата система с Н, а определено състояние на системата с n-мерния вектор X = (X₁, X₂, ... X_n), където n е броя на степените на свобода. Множеството от възможни състояния ={X} образува т.нар. фазово пространство. Тогава наблюдаемата величина А ще се определи от всички възможни състояния на системата с помощта на т.нар. функция на разпределение f(H(X)):

К
ъдето Z е т.нар. статистически интеграл, т.е. А се получава чрез усредняване по статистически ансамбъл (по състояния). В КЕ много по-често се усреднява по траектория във фазовото пространство т.е. по време.

Хипотезата за ергодичност допуска, че <A>=A_. Така усредняването по време за безкраен интервал се равнява на усредняването по състояния тъй като се предполага че за безкрайно време системата ще премине през всички възможни състояния. На практика очевидно е невъзможно времето да е безкрайно затова хипотезата за ергодичност е приблизително вярна за крайно време: <A>A_t

Има два основни метода за преместване на системата във фазовото пространство:

1) детерминистичен

2) вероятностен (стохастичен)

Идеята на детерминистичния метод е, че се използва собствената динамика на системата т.е. задават се уравненията на движение и се интегрират във времето. Наборът от частици на системата се описва с помощта на класическата механика т.е. при зададени начални координати X^N(0) = {X₁(0), X₂(0), ... X_N(0)} и моменти p^N(0) = {p₁(0), p₂(0), ... p_N(0)} трябва да намерим траекторията във фазовото пространство (X^N(t), p^N(t)). Най-популярният представител на детерминистичния подход е метода за компютърно симулиране наречен Молекулярна Динамика (МД).

Методът МД е въведен от Алдер и Вайнрайт през 1950 г. за изучаване на прости течности чрез модел на взаимодействие между твърди сфери. Изследван е фазовия преход течност-твърда фаза.

Следващият успех е през 1964 г. когато Раман прави симулация на течен аргон като използва реалистичен потенциал на взаимодействие (Ленард-Джоунс) в система от 864 атома. През 1974 се прави първата симулация на течна вода.

Първата симулация на протеин е през 1977 г. Днес всички белтъчни биополимери се моделират с помощта на МД. Други приложения са: рентгенова кристалография, неутронен магнитен резонанс и др. за определяне структурата на материалите.

При стохастичните методи подходът е друг. Ясно е, че трябва да се изчисляват различните конфигурации (състояния) на системата. Проблемът е как да придвижим системата от едно състояние в друго. Това става чрез вероятностна еволюция, която се осъществява чрез марковски процес, който е вероятностен аналог на собствената динамика при детерминистичния метод. Затова стохастичният метод може да моделира системи без собствена динамика. Най-популярният представител на стохастичния подход е метода Монте Карло (МК), при който симулацията се извършва с помощта на поредица от случайни числа.

За разлика от МД методите МК не се нуждаят от дефиниране на диференциални уравнения, описващи поведението на системата. Единственото изискване е физическата система да бъде описана с помощта на една или повече функции на разпределение на вероятността (ФРВ). След като знаем набора от ФРВ за дадена система можем да пресъздадем множество варианти на еволюция на системата с помощта на метода МК и след усредняване по състоянията да получим резултати за различни макровеличини описващи системата. Предимство на метода е, че може да се предскаже систематичната грешка и оттам да се оцени необходимия брой наблюдения за постигане на дадена грешка.

Приложенията на метода МК са многобройни и разнородни и обхващат широк кръг от проблеми, Някои от тези приложения са симулации на:

• 
  ядрени реакции
  
• 
  квантова хромодинамика
  
• 
  терапия на рак чрез облъчване
  
• 
  транспортен поток
  
• 
  звездни еволюции
  
• 
  предсказване на индекса Доу-Джонс
  
• 
  изследване на нефтени находища
  
• 
  проектиране на интегрални схеми и т.н.
  

1. 
   Първото ограничение е свързано с крайното време за
   

2. 
   Второто основно ограничение е свързано с крайния размер на
   

3. 
   Алгоритмите за числено интегриране не са безкрайно точни:
   

4. 
   В методите МД междумолекулярните сили не са точно известни.
   

5. 
   В методите МК основното ограничение е свързано с качествата на
   

Компютърни алгоритми

Преди да може да се реши даден проблем на компютър трябва да се инструктира компютъра какви конкретни операции да извърши. Процесът включва две стъпки: Първо, трябва да се преобразува проблема, зададен обикновенно като уравнение, в набор от логически стъпки, които компютърът може да изпълни. Това е т.нар. числов компютърен алгоритъм. Второ, трябва да се информира компютъра как да проследи и изпълни тези стъпки. Това се реализира с помощта на някакъв компютърен език.

Щ
е разгледаме алгоритъма на една проста програма, описваща 1-D движение на една частица с маса m под действието на хармонична сила. Нека опишем движението на частицата с уравнението на Нютон: F=ma, където а е ускорението. Ако разделим времето на малки равни интервали  можем да определим средната скорост и средното ускорение за интервала [t_n , t_n+1 ]:

Най-простият алгоритъм за определяне на положението и скороста на частицата в момент t_n+1 от известните стойности в предишния момент t_n се получава след преобразуване на горните уравнения:

В
най-дясната част на последното уравнение сме приели m=1, а за хармоничната външна сила се изпълнява условието f(x)=-x. Очевидно ако знаем началните условия за скороста и координатата на частицата можем да определим тяхните стойности във всеки по-късен момент с помощта на показания рекурсивен алгоритъм наричан метод на Ойлер.

Този прост пример показва как се конструират повечето алгоритми. Първо, физическите уравнения се привеждат в дискретна форма т.е. в уравнения на разликите. След това желаните величини или решения се изразяват по рекурсивен начин т.е. стойностите в по-късен момент се изразяват чрез стойностите в предишен момент.

Логическите стъпки в един алгоритъм могат да бъдат последователни, паралелни или итеративни. Как да се използват свойствата на дадена задача при конструирането на бърз и точен алгоритъм е много важен момент в изчислителната наука. Вторият важен момент е да се напише чиста и добре организирана компютърна програма.

Зад. 1.1: Заредете файл с програмата написана с език С, фортран или JAVA. Kомпилирайте и стартирайте програмата, написана според горния алгоритъм. Резултатите за x_n, v_n и t_n се извеждат на екрана след всеки MMAX итерации.

• 
  Модифицирайте програмата така, че да се въвеждат началните условия от клавиатурата и да се записват резултатите в три текстови файла за x, v и t.
  
• 
  Постройте графични зависимости x(t) и v(t) като използвате графичните възможности на друга програма, примерно Excel, MatLab или Origin.
  

Една класическа числова задача е сумирането на ред за да се оцени функцията, която е разложена в ред. Ще разгледаме като пример степенния ред на експоненциалната функция.

Искаме да използваме горния ред за да изчислим exp(x) за х=0.1, 1, 10, 100 и 1000 с абсолютна грешка не надвишаваща 10^-8. Как обаче да разберем кога да спрем сумирането? Лесно е да допуснем че грешката в сумирането се определя от последния член в реда. Тогава можем да приемем правилото че сумирането спира когато се изпълни условието:

където term е последния член в реда, а sum е натрупаната сума.

• 
  Нека се запише във файл таблица на стойностите на
  

• 
  Проверете програмата при големи стойности на х, примерно 50. Оценете грешката, като приемем че вградената функция дава точен резултат. Наблюдава се голяма грешка дължаща се на вадене на две големи, но близки по стойност числа.
  
• 
  Намалява ли грешката, ако пресметнем exp(x) и след това вземем реципрочната стойност 1/exp(x). В този случай нямаме вадене на големи числа, а само събиране.
  

// An example of studying the motion of a particle in

// one dimension under an elastic force.
import java.lang.*;

public class Motion {

static final int n = 100000, j = 500;

public static void main(String argv[]) {

double x[] = new double[n+1];

double v[] = new double[n+1];

double dt = 2*Math.PI/n;

x[0] = 0;

v[0] = 1;

for (int i=0; i

x[i+1] = x[i]+v[i]*dt;

v[i+1] = v[i]-x[i]*dt;

}
// Output the result in every j time steps

double t = 0;

double jdt = j*dt;

for (int i=0; i<=n; i+=j) {

System.out.println(t +" " + x[i] + " " + v[i]);

t += jdt;

}

}

}

term=1, sum=1, eps=10^(-8)

do

term=-term*x/i

while abs(term/sum)>eps