Форум программистов
 

Восстановите пароль или Зарегистрируйтесь на форуме, о проблемах и с заказом рекламы пишите сюда - alarforum@yandex.ru, проверяйте папку спам!

Вернуться   Форум программистов > IT форум > Помощь студентам
Регистрация

Восстановить пароль
Повторная активизация e-mail

Купить рекламу на форуме - 42 тыс руб за месяц

Ответ
 
Опции темы Поиск в этой теме
Старый 30.07.2010, 23:45   #1
FedyaninT
Новичок
Джуниор
 
Регистрация: 30.07.2010
Сообщений: 2
По умолчанию Помогите с Mathlab'ом

Доброго времени суток.
Прошу помощи у знающих людей. Сам я не знаком с матлабом, но тут очнеь понадобилось разобраться с кодом небольшой программы (нужно реализовать ее на C). Пожалуйсто прокоментируйте некоторые моменты (выделил жирным), что там происходит (ну такие вещи как x.^2 вроде понятно, возведение в степень, а вот x^*y уже вызывают затруднения), а то синтаксис языка совсем не интуитивно понятный.

Код:
% Spatial grid and initial condition:
N = 128;
x = 32*pi*(1:N)’/N;
u = cos(x/16).*(1+sin(x/16));
v = fft(u);
% Precompute various ETDRK4 scalar quantities:
h = 1/4; % time step
k = [0:N/2-1 0 -N/2+1:-1]’/16; % wave numbers
L = k.^2 - k.^4; % Fourier multipliers
E = exp(h*L); E2 = exp(h*L/2);
M = 16; % no. of points for complex means
r = exp(1i*pi*((1:M)-.5)/M); % roots of unity
LR = h*L(:,ones(M,1)) + r(ones(N,1),:);
Q = h*real(mean( (exp(LR/2)-1)./LR ,2));
f1 = h*real(mean( (-4-LR+exp(LR).*(4-3*LR+LR.^2))./LR.^3 ,2));
f2 = h*real(mean( (2+LR+exp(LR).*(-2+LR))./LR.^3 ,2));
f3 = h*real(mean( (-4-3*LR-LR.^2+exp(LR).*(4-LR))./LR.^3 ,2));
% Main time-stepping loop:
uu=u;tt=0;
tmax = 150; nmax = round(tmax/h); nplt = floor((tmax/100)/h);
g = -0.5i*k;
for n = 1:nmax
t = n*h;
Nv = g.*fft(real(ifft(v)).^2);
a = E2.*v + Q.*Nv;
Na = g.*fft(real(ifft(a)).^2);
b = E2.*v + Q.*Na;
Nb = g.*fft(real(ifft(b)).^2);
c = E2.*a + Q.*(2*Nb-Nv);
Nc = g.*fft(real(ifft(c)).^2);
v = E.*v + Nv.*f1 + 2*(Na+Nb).*f2 + Nc.*f3;
if mod(n,nplt)==0
u = real(ifft(v));
uu = [uu,u]; tt = [tt,t];
end
end
% Plot results:
surf(tt,x,uu), shading interp, lighting phong, axis tight
view([-90 90]), colormap(autumn); set(gca,’zlim’,[-5 50])
light(’color’,[1 1 0],’position’,[-1,2,2])
material([0.30 0.60 0.60 40.00 1.00]);
Заранее благодарен!
FedyaninT вне форума Ответить с цитированием
Старый 02.08.2010, 01:20   #2
arcer
Пользователь
 
Регистрация: 26.01.2010
Сообщений: 42
По умолчанию

Цитата:
Сообщение от FedyaninT Посмотреть сообщение
Доброго времени суток.
Прошу помощи у знающих людей. Сам я не знаком с матлабом, но тут очнеь понадобилось разобраться с кодом небольшой программы (нужно реализовать ее на C). Пожалуйсто прокоментируйте некоторые моменты (выделил жирным), что там происходит (ну такие вещи как x.^2 вроде понятно, возведение в степень, а вот x^*y уже вызывают затруднения), а то синтаксис языка совсем не интуитивно понятный.

Код:
% Spatial grid and initial condition:
N = 128;
x = 32*pi*(1:N)’/N;
u = cos(x/16).*(1+sin(x/16));
v = fft(u);
% Precompute various ETDRK4 scalar quantities:
h = 1/4; % time step
k = [0:N/2-1 0 -N/2+1:-1]’/16; % wave numbers
L = k.^2 - k.^4; % Fourier multipliers
E = exp(h*L); E2 = exp(h*L/2);
M = 16; % no. of points for complex means
r = exp(1i*pi*((1:M)-.5)/M); % roots of unity
LR = h*L(:,ones(M,1)) + r(ones(N,1),:);
Q = h*real(mean( (exp(LR/2)-1)./LR ,2));
f1 = h*real(mean( (-4-LR+exp(LR).*(4-3*LR+LR.^2))./LR.^3 ,2));
f2 = h*real(mean( (2+LR+exp(LR).*(-2+LR))./LR.^3 ,2));
f3 = h*real(mean( (-4-3*LR-LR.^2+exp(LR).*(4-LR))./LR.^3 ,2));
% Main time-stepping loop:
uu=u;tt=0;
tmax = 150; nmax = round(tmax/h); nplt = floor((tmax/100)/h);
g = -0.5i*k;
for n = 1:nmax
t = n*h;
Nv = g.*fft(real(ifft(v)).^2);
a = E2.*v + Q.*Nv;
Na = g.*fft(real(ifft(a)).^2);
b = E2.*v + Q.*Na;
Nb = g.*fft(real(ifft(b)).^2);
c = E2.*a + Q.*(2*Nb-Nv);
Nc = g.*fft(real(ifft(c)).^2);
v = E.*v + Nv.*f1 + 2*(Na+Nb).*f2 + Nc.*f3;
if mod(n,nplt)==0
u = real(ifft(v));
uu = [uu,u]; tt = [tt,t];
end
end
% Plot results:
surf(tt,x,uu), shading interp, lighting phong, axis tight
view([-90 90]), colormap(autumn); set(gca,’zlim’,[-5 50])
light(’color’,[1 1 0],’position’,[-1,2,2])
material([0.30 0.60 0.60 40.00 1.00]);
Заранее благодарен!
Значит так, без претензий на 100% достоверность (давно работал, забыл немного):
x = 32*pi*(1:N)’/N; - (1:N) - N - мерный вектор, заполненый 1(погугли), (1:N)’ - столбик, выражение, думаю, понятно.
u = cos(x/16).*(1+sin(x/16)); - .* - операция поэлементного умножения векторов (т.е. [15 2].*[2 13] == [20 26]), не забываем, что х - вектор
[0:N/2-1 0 -N/2+1:-1] - погугли
nplt = floor((tmax/100)/h); - округление, наибольшее целое, которое меньше или равно аргументу
mod(n,nplt) - взятия остатка от деления n/nplt
u = real(ifft(v)); - гуглофон
uu = [uu,u]; - конкатенация, кажеться

все остальное - одно из вышеупомянутых.
arcer вне форума Ответить с цитированием
Старый 02.08.2010, 19:31   #3
FedyaninT
Новичок
Джуниор
 
Регистрация: 30.07.2010
Сообщений: 2
Радость

Большое спасибо arcer.
FedyaninT вне форума Ответить с цитированием
Ответ


Купить рекламу на форуме - 42 тыс руб за месяц