TimeCoder2

Здесь будут размещаться программы/скрипты на разных языках, реализующие PaPuRi-алгоритм. Складывается впечатление что исследователи обсуждающие конкретные алгоритмы математической физики (не в превдокодах, а в виде, пригодном для трансляции и исполнения на вычислительной технике) очень часто используют для этого MATLAB. Однако, тут упоминаются и бесплатные альтернативы програм/сред трансляции подобных алгоритмов.

Что делает любой из этих кодов в самом простом случае, какую задачу решает?
(есть и обобщения, но начнём с простого) Надо численно решить конкретную систему из 2 линейных дифференциальных уравнений в частных производных с заданными начальными и граничными условиями.
di/dt+du/dx+Ri=0; du/dt+di/dx+Gu=0,
где u на концах равна нулю и обе величины i,u заданы в начальный момент времени. Простейший случай — u на внутреннем интервале равно единице, остальные начальные данные равны нулю. При этом, даже для R=G=0, при известном науке точном, аналитическом решении, известный численный метод FDTD ведёт себя плохо. Мы не будем распространяться тут о смысле данной задачи для молниеотводов и линий электропередач, или даже о выводе алгоритма; уже на одном этом сайте в разделе PaPuRi Group вы найдёте об этом массу публикаций и книг.

Итак, реализации алгоритма.

MATLAB (не векторизованный):
tic; clear; n=200; nn=n+1; h=1; tau=h; k=400;
for i=1:n X(i)=i*h-h/2; R(i)=0; G(i)=0; D(i)=0; U(i)=0; end;
for i=1:n
al(i)=(R(i)+G(i))/2; AD(i)=1+tau*al(i);
AD1(i)=1+tau*(al(i)-R(i));
AU1(i)=1+tau*(al(i)-G(i)); end;
for i=76:125 U(i)=1; end;
UA(1)=0;UA(nn)=0;
for j=1:k
for i=2:n
DA(i)=(D(i-1)+D(i)+U(i-1)-U(i))/2;
UA(i)=(U(i-1)+U(i)+D(i-1)-D(i))/2; end;
DA(1)=D(1)-U(1); DA(nn)=D(n)+U(n);
for i=1:n
D(i)=(UA(i)-UA(i+1)+D(i)*AD1(i))/AD(i);
U(i)=(DA(i)-DA(i+1)+U(i)*AU1(i))/AD(i); end;
plot(X,D,’k’,X,U,’r’); pause(0.001);
end; toc;

Проще понять тем кто не знаком с векторными операциями. На Octave (бесплатном аналоге MATLAB), работает слишком медленно, тем более что задержку приходится увеличивать до 0.1 для того чтобы он вообще успевал что-то отрисовывать.

MATLAB (векторизованный):
tic; clear; n=200; nn=n+1; h=1; tau=h; k=400;
X=h/2:h:n*h-h/2; R=zeros(1,n); G=zeros(1,n); D=zeros(1,n); U=zeros(1,n); U(76:125)=1;
al=(R+G)/2; AD=1+tau*al; AD1=1+tau*(al-R); AU1=1+tau*(al-G);
for j=1:k
DA=[D(1)-U(1) (D(1:n-1)+D(2:n)+U(1:n-1)-U(2:n))/2 D(n)+U(n)];
UA=[0 (U(1:n-1)+U(2:n)+D(1:n-1)-D(2:n))/2 0];
D=(UA(1:n)-UA(2:nn)+D.*AD1)./AD;
U=(DA(1:n)-DA(2:n+1)+U.*AU1)./AD;
plot(X,D,’k’,X,U,’r’); pause(0.001);
end; toc;

Компактнее, в силу чего проще понять знакомым с векторными операциями (а иначе не зря ли учили MATLAB?) На Octave вычисления идут куда быстрее. Странно что из-за операций извлечения векторов без первого либо последнего элемента, сам MATLAB работает чуть медленнее, что, впрочем, не особо заметно на фоне операций рисования и задержки.

MATLAB (с пропуском кадров):
tic; clear; n=200; nn=n+1; h=1; tau=h; k=400;
X=h/2:h:n*h-h/2; R=zeros(1,n); G=zeros(1,n); D=zeros(1,n); U=zeros(1,n); U(76:125)=1;
al=(R+G)/2; AD=1+tau*al; AD1=1+tau*(al-R); AU1=1+tau*(al-G);
for j=1:k
DA=[D(1)-U(1) (D(1:n-1)+D(2:n)+U(1:n-1)-U(2:n))/2 D(n)+U(n)];
UA=[0 (U(1:n-1)+U(2:n)+D(1:n-1)-D(2:n))/2 0];
D=(UA(1:n)-UA(2:nn)+D.*AD1)./AD;
U=(DA(1:n)-DA(2:n+1)+U.*AU1)./AD;
if(mod(j,10)==0) disp(sprintf(‘%d’,j)); plot(X,D,’k’,X,U,’r’); pause(0.1); end;
end; toc;

При большом количестве итераций бывает скучно смотреть на каждый кадр, поэтому будем показывать, например, каждый 10й.

MATLAB (с вызовом модуля написанного на C):
tic; clear; n=200; h=1; tau=h; k=400; period=1; delay=0.001; j=0;
X=h/2:h:n*h-h/2; R=zeros(1,n); G=zeros(1,n); D=zeros(1,n); U=zeros(1,n); U(76:125)=1;
xfrm(D,U,R,G,tau,k,period,j); toc;

А вот этот код работает во много раз быстрее, благодаря тому что главная работа переписана на C. Главное было переписать хотя бы вычислительную часть внутри цикла по j, но в итоге оказалось целесообразнее вызывать команду рисования уже из C, переписав таким образом все вычисления. Кроме инициализации, которая может быть удобнее для пользователей на MATLAB, и по скорости не критична. Для компиляции использовалась команда mex xfrm.c

xfrm.c: модуль на C для вызова из MATLAB
#include «mex.h»
void mexFunction(int nlhs,mxArray *plhs[],int nrhs,const mxArray *prhs[])
{int n,i,j,k,s,ss,period,delay;
double al,tau;
double* D,* U,* R,* G,* DA,* UA,* AD1,* AU1,* AD;
n=mxGetN(prhs[0]); D=mxGetPr(prhs[0]); U=mxGetPr(prhs[1]);
R=mxGetPr(prhs[2]); G=mxGetPr(prhs[3]);
tau=mxGetPr(prhs[4])[0]; k=mxGetPr(prhs[5])[0]; period=mxGetPr(prhs[6])[0];
s=n*sizeof(double); ss=(n+1)*sizeof(double);
AD1=malloc(s); AU1=malloc(s); AD=malloc(s); DA=malloc(ss); UA=malloc(ss);
for(i=0;i<n;i++)
{al=(R[i]+G[i])/2;
AD[i]=1+tau*al;
AD1[i]=1+tau*(al-R[i]);
AU1[i]=1+tau*(al-G[i]);
}
UA[0]=0; UA[n]=0;
for(j=0;j<k;j++)
{for(i=1;i<n;i++)
{DA[i]=(D[i-1]+D[i]+U[i-1]-U[i])/2;
UA[i]=(U[i-1]+U[i]+D[i-1]-D[i])/2;
}
DA[0]=D[0]-U[0]; DA[n]=D[n-1]+U[n-1];
for(i=0;i<n;i++)
{D[i]=(UA[i]-UA[i+1]+D[i]*AD1[i])/AD[i];
U[i]=(DA[i]-DA[i+1]+U[i]*AU1[i])/AD[i];
}
if((j+1)%period==0)
{mxGetPr(prhs[7])[0]=j+1;
mexEvalString(«plot(X,D,’k’,X,U,’r’);title([‘j=’,num2str(j)]);pause(delay);»);
}
}
free(UA); free(DA); free(AD); free(AU1); free(AD1);
}

Последнее навевает вопрос — а зачем вообще тогда нужен MATLAB для данной задачи, кроме простой визуализации (в языках типа C нет для этого универсального встроенного инструмента), и, возможно, простоты инициализации?

GNUPLOT:
n=200; nn=n+1; h=1; tau=h; k=400;
array X[n]; array R[n]; array G[n]; array D[n]; array U[n];
array AD[n]; array AD1[n]; array AU1[n]; array DA[nn]; array UA[nn];
do for [i=1:n]{
X[i]=i*h-h/2.0; R[i]=0; G[i]=0; D[i]=0; U[i]=0;
al=(R[i]+G[i])/2.0; AD[i]=1+tau*al; AD1[i]=1+tau*(al-R[i]); AU1[i]=1+tau*(al-G[i]);
}
do for [i=76:125]{U[i]=1};
set linetype 1 lc «red» lw 8; set linetype 2 lc «black» lw 5;
UA[1]=0; UA[nn]=0;
do for [j=1:k]{
do for [i=2:n]{DA[i]=(D[i-1]+D[i]+U[i-1]-U[i])/2.0; UA[i]=(U[i-1]+U[i]+D[i-1]-D[i])/2.0;}
DA[1]=D[1]-U[1]; DA[nn]=D[n]+U[n];
do for [i=1:n]{D[i]=(UA[i]-UA[i+1]+D[i]*AD1[i])/AD[i]; U[i]=(DA[i]-DA[i+1]+U[i]*AU1[i])/AD[i];}
if(j%1==0){set title sprintf(‘Step %d’,j); plot U with lines,D with lines; pause 0.001;}
}

Решение было найдено следующее: выяснилось что бесплатная программа GnuPlot содержит в себе вполне достаточный инструментарий не только для рисования графиков, но и для проведения соответствующих расчётов; работает достаточно эффективно.

JavaScript:

TimeCoder — А на этой странице вы можете ознакомиться с реализацией алгоритма вообще без необходимости инсталлировать какой-либо конкретный специализированный софт 🙂 Вычисления идут прямо в браузере! К сожалению, с инициализацией данных в этом последнем случае сложнее, но приведены несколько примеров, и даже начинающему программисту не должно составить труда сгенерировать надлежащий текст из чисел, описывающий интересующую его задачу, на своём любимом языке программирования, а задачи крошечных размеров для ознакомления и вовсе можно ввести вручную.

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *