Материал: Matlab. Практический подход. Самоучитель

Самоучитель Matlab

function s = u2(A,x,y,t,k) s=0;

Nmax=length(A); for m=1:Nmax

for n=1:Nmax s=s+A(m,n)*cos(pi*t*sqrt(m^2+n^2*k^2)).*sin(pi*m*x).*sin(pi*n*y); end

end end

Функцией u2() вычисляется ряд для решения дифференциального уравнения. Особенность функции, по сравнению с предыдущим случаем, состоит в том, что, помимо пространственных переменных x и y, времени t и параметра k, первым аргументом функции передается список A коэффициентов разложения функции начального профиля по базисным функциям. Количество слагаемых в ряде определяется количеством элементов (коэффициентов) в списке A. Коэффициенты этого списка вычисляются с помощью другой функции, которая называется Coefs(). Ее код приведен ниже:

function A=Coefs(f) Nmax=30; A=zeros(Nmax,Nmax); for m=1:Nmax

for n=1:Nmax F=@(x,y)(f(x,y).*sin(pi*m*x).*sin(pi*n*y)); A(m,n)=4*quad2d(F,0,1,0,1);

end

end end

Аргумент у функции один – указатель f на функцию, определяющую профиль начального распределения мембраны. Предполагается, что это функция двух переменных. В теле функции Coefs() определяется переменная Nmax. Ее значение в данном случае равно 100, и это количество слагаемых ряда (по каждой из индексных переменных). Командой A=zeros(Nmax,Nmax) определяется нулевая квадратная матрица, в которую будут заноситься вычисленные коэффициенты разложения. В теле вложенных операторов цикла на основе указателя f командой F=@(x,y)(f(x,y).*sin(pi*m*x).*sin(pi*n*y))

определяется подынтегральная функция, на основе которой вычисляются коэффициенты ряда. Коэффициенты вычисляются командой A(m,n)=4*quad2d(F,0,1,0,1). В этой команде использована функция quad2d(), предназначенная для вычисления интегралов по плоской области. Первым аргументом функции передается указатель на подын-

306

Глава 7. Уравнения математической физики

тегральную функцию. Другие аргументы – пределы интегрирования по каждой из переменных интегрирования.

Для решения задачи создаем файл Membrane2.m со следующим кодом:

%Функция начального профиля: f=@(x,y)exp(-200*((x-0.5).^2+(y-0.5).^2));

%Вычисление коэффициентов:

A=Coefs(f);

%Количество кадров (+1): N=300;

%Продолжительность по времени: T=3;

%Пространственные координаты: x=0:0.02:1;

y=0:0.02:1;

%Геометрический фактор:

k=1;

%Матрицы для графических утилит: [X,Y]=meshgrid(x,y);

for i=0:N

%Момент времени:

t=i*T/N;

%Расчет профиля мембраны: Z=u2(A,X,Y,t,k);

%Поверхность:

surf(X,Y,Z);

%Диапазон значений по осям: axis([0 1,0 1,-0.4 0.4]);

%Координатная сетка:

grid on;

%Кадр анимации: F(i+1)=getframe; end

%Отображение анимации: movie(F);

Функцию начального профиля создаем командой f=@(x,y)exp(-200* ((x-0.5).^2+(y-0.5).^2)). Здесь речь идет о функциональной зависимости f (x,y) = exp(−200(x2 + y2)). Это локализованный в центре квадратной мембраны пик Гауссового типа. В начальный момент мембрана (ее профиль) имеет вид, как на рис. 7.32.

На заметку

Желающие могут подумать, как с минимальными затратами времени и ресурсов построить такой график.

307

Самоучитель Matlab

Рис. 7.32. Начальный профиль мембраны с пиком в центре

Некоторые кадры анимации, созданной в результате выполнения приведенного выше кода, показаны в табл. 7.2.

Табл. 7.2. Колебания мембраны в случае начального профиля в виде пика

308

Глава 7. Уравнения математической физики

309

Самоучитель Matlab

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

На заметку

Выше был использован подход, при котором коэффициенты разложения ряда вычислялись вне тела функции-решения уравнения. Это было сделано не случайно. Дело в том, что если ряд вычисляется по большому количеству слагаемых, значительное время уходит на расчет этих слагаемых. Если коэффициенты вычислять в теле функции, то такие расчеты будут выполняться каждый раз, когда вызывается функция решения. В использованном выше подходе коэффициенты были вычислены единожды и записаны в матрицу. А при вызове функции решения коэффициенты считываются из этой матрицы. С другой стороны, если матрица очень большая, могут возникнуть проблемы с системными ресурсами, в то время как расчет коэффициентов непосредственно при вызове функции решения можно реализовать так, что использование системных ресурсов сведено к минимуму (зато это будет долго по времени). Тот или иной подход обычно выбирают, исходя из конкретики решаемой задачи и наличия ресурсов для ее решения. В любом случае хорошо, что есть из чего выбирать.

310