Как уменьшить время счета

@Golem_ru · Регистрация: 20.06.2013

Author24 — интернет-сервис помощи студентам

Ребят, есть программа, проблема в том, что создавая в цикле анонимную функцию, время счета программы очень велико.
71 строчка -Создание анонимной функции F1=@ (x) 1./sqrt(1-k^2.*(1-sin(x).^2)), k вычисляется в цикле.
72 строчка - Вычисление интеграла K=quad(F1,0,pi/2);
Для проверки счета количество N вводить 10, больше не надо.
Считал интеграл и методом Симпсона, но время счета было ещё больше.

Помогите разобраться с этим циклом, чтобы время счета было гораздо меньше. Спасибо заранее за помощь

Matlab M

clc
clear all
close all
e0 = 8.85e-12; %[Ф/м]
N=input('Введите количество эквивалентных зарядов на 1/4 окружности фланца N= ');
 
R_f=0.05;  %cm
R0=5; %cm
U=1; %кВ
deltamax=1e-6;
e=4; % УТОЧНИТЬ!!!!!!!!!
R_t=0.5;
tic
delta_alfa=pi/N;
L=delta_alfa*R_f/(1+delta_alfa*2);
H=2*L;
L_f=2*R0; %Меняем
L_t=L_f;
l=ceil((L_f-R_f)/L);
Y(1:2*l+N+2)=0;
Y(2*l+N+3:3*l+N+4)=U;
 
ri=zeros(1,3*l+N+4);
zi=zeros(1,3*l+N+4);
rki=zeros(1,3*l+N+4);
zki=zeros(1,3*l+N+4);
%Координаты ЭЗ НА фланце 
for i=N/2+1:N;
    ri(i)=R0+R_f+(R_f-H)*cos(delta_alfa*(i));
    zi(i)=(R_f-H)*sin(delta_alfa*(i));
    rki(i)=R0+R_f+R_f*cos(delta_alfa*(i));
    zki(i)=R_f*sin(delta_alfa*(i));
end
zi(1:N/2)=-zi(N/2+1:N);
ri(1:N/2)=ri(N/2+1:N);
rki(1:N/2)=rki(N/2+1:N);
zki(1:N/2)=-zki(N/2+1:N); 
 
for i=N+1:l+N+1;
    ri(i)=R0+R_f+L*(i-N-1);   
end
zi(N+1:l+N+1)=R_f-H;
rki(N+1:l+N+1)=ri(N+1:l+N+1);
zki(N+1:l+N+1)=R_f;
 
ri(l+N+2:2*l+N+2)=ri(N+1:l+N+1) ;
zi(l+N+2:2*l+N+2)=-zi(N+1:l+N+1);
rki(l+N+2:2*l+N+2)=rki(N+1:l+N+1);
zki(l+N+2:2*l+N+2)=-zki(N+1:l+N+1);
%Координаты ЭЗ НА трубе
 
for i=2*l+N+3:3*l+N+4;
    zi(i)=-L_t/2+L*(i-2*l-N-3);
end
ri(2*l+N+3:3*l+N+4)=R_t-H;
zki(2*l+N+3:3*l+N+4)=zi(2*l+N+3:3*l+N+4);
rki(2*l+N+3:3*l+N+4)=R_t;
 
tic
 
 
% A=zeros(3*l+N+4,3*l+N+4);
B=zeros(1,3*l+N+4);
 for I=1:3*l+N+4;
     for J=1:3*l+N+4;
         a=1/(sqrt((zki(J)-zi(I))^2+(rki(J)+ri(I))^2));
         k=sqrt((4*rki(J)*ri(I))/((zki(J)-zi(I))^2+(rki(J)+ri(I))^2));
 
 
 
%          F1=@ (x) 1./sqrt(1-k^2.*(1-sin(x).^2));  
%          K=quad(F1,0,pi/2);
%          A(I,J)=a*K;
     end
       B(I)=Y(I)*(2*pi^2*e0*e);
 end
  t=toc
%  C=B';
%  T=A^(-1)*C;

Добавлено через 17 минут
Поправка 71 строчки

Matlab M

1	% F1=@ (x) 1./sqrt(1-k^2.*sin(x).^2);

Зосима · 08.10.2015, 21:19

А как насчет численного интегрирования? оно всяко быстрее будет

Matlab M

B=zeros(1,3*l+N+4);
x = 0:0.01:pi/2;
for I=1:3*l+N+4;
    for J=1:3*l+N+4;
        a=1/(sqrt((zki(J)-zi(I))^2+(rki(J)+ri(I))^2));
        k=sqrt((4*rki(J)*ri(I))/((zki(J)-zi(I))^2+(rki(J)+ri(I))^2));
        F1 = 1./sqrt( 1 - k^2.*sin(x).^2 );
        K = trapz(x,F1);
        A(I,J) = a*K;
    end
    B(I)=Y(I)*(2*pi^2*e0*e);
end

@Golem_ru · 08.10.2015, 22:32 **[ТС]**

Спасибо, сам допер до этого, только вопрос, достаточно точный расчет то получится методом трапеций?? Или лучше сделать метод Симпсона, просто реализовывая метод Симпсона, фигурирует опять анонимная функция, что не очень приятно по времени счета.

Зосима · 09.10.2015, 11:48

Golem_ru, не, для 9млн. элементов (при N=10) метод Симпсона мягко говоря не очень хорошая идея

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

(полностью не выйдет, т.к. памяти не хватит)
Но чтобы размеры векторов x и a,k были равны нужно использовать meshgrid, т.е. из двух векторов разной длинны сделать матрицы одинакового размера:

Matlab M

A = zeros(3*l+N+4);
B=zeros(1,3*l+N+4);
x = 0:0.01:pi/2;
for I=1:3*l+N+4;
    [X,a] = meshgrid(x, 1./(sqrt((zki-zi(I)).^2+(rki+ri(I)).^2)) );
    [X,k] = meshgrid(x, sqrt((4*rki*ri(I))./((zki-zi(I)).^2+(rki+ri(I)).^2)) );
    A(I,:) = trapz(x, 1./sqrt(1-k.^2.*sin(X).^2), 2); % интегрируем по столбцам
    B(I)=Y(I)*(2*pi^2*e0*e);
end

для одной итерации (for I=1 ...) время на моем тракторе составило t=0.25c, т.е. для обработки 3200 элементов (длинна массива I) понадобится ~800c ~ 1мин 20с, хотя спинным мозгом чувствую, что на деле будет больше

При желании можно прикрутить сюда waitbar

Matlab M

h = waitbar(0,'Идет рассчет...');
for I=1:3*l+N+4;
    [X,a] = meshgrid(x, 1./(sqrt((zki-zi(I)).^2+(rki+ri(I)).^2)) );
    [X,k] = meshgrid(x, sqrt((4*rki*ri(I))./((zki-zi(I)).^2+(rki+ri(I)).^2)) );
    A(I,:) = trapz(x, 1./sqrt(1-k.^2.*sin(X).^2), 2); % интегрируем по столбцам
    B(I)=Y(I)*(2*pi^2*e0*e);
    waitbar( I/(3*l+N+4) )
end
close(h)

*точность метода трапеций можно оценить с помощью правила Рунге, т.е. посчитать интеграл для одной конкретной пары значений a, k с разными шагами переменной х:

Matlab M

a = 1/(sqrt((zki(1)-zi(1))^2+(rki(1)+ri(1))^2));
k = sqrt((4*rki(1)*ri(1))/((zki(1)-zi(1))^2+(rki(1)+ri(1))^2));
 
x = 0:0.01:pi/2; % считаем с шагом 0.01
F1 = trapz(x, 1./sqrt(1 - k^2.*sin(x).^2) );
x = 0:0.005:pi/2; % считаем с шагом 0.005
F2 = trapz(x, 1./sqrt(1 - k^2.*sin(x).^2) );
R = 1/3*abs(F2-F1) % абс. точность при шаге 0.01

@Nick07 · 11.10.2015, 22:34

1. Поставить опер. памяти по максимуму.
2. Пробовать parfor:
http://www.mathworks.com/help/... arfor.html
3. Пробовать spmd:
http://www.mathworks.com/help/distcomp/spmd.html
4. Пробовать GPU:
http://www.mathworks.com/help/... array.html

Добавлено через 36 минут
http://matlab.exponenta.ru/for... 19631.html

@Golem_ru · 12.10.2015, 12:49 **[ТС]**

Про waitbar не знал, удобная штука))) Да, ребят такой вопрос, как выделить всю оперативную память матлаб? Много где читал что по умолчанию 512 мб стоит, а когда идут тяжелые вычисления он мол всю задействует, но смотря как он считает, я так думаю что 512 и стоит.

@Golem_ru 6 / 4 / 1 Регистрация: 20.06.2013 Сообщений: 24
	08.10.2015, 22:32 [ТС]	3
	Спасибо, сам допер до этого, только вопрос, достаточно точный расчет то получится методом трапеций?? Или лучше сделать метод Симпсона, просто реализовывая метод Симпсона, фигурирует опять анонимная функция, что не очень приятно по времени счета. 0

@Nick07 533 / 438 / 47 Регистрация: 17.07.2013 Сообщений: 2,236
	11.10.2015, 22:34	5
	1. Поставить опер. памяти по максимуму. 2. Пробовать parfor: http://www.mathworks.com/help/... arfor.html 3. Пробовать spmd: http://www.mathworks.com/help/distcomp/spmd.html 4. Пробовать GPU: http://www.mathworks.com/help/... array.html Добавлено через 36 минут http://matlab.exponenta.ru/for... 19631.html 1

@Golem_ru 6 / 4 / 1 Регистрация: 20.06.2013 Сообщений: 24
	12.10.2015, 12:49 [ТС]	6
	Про waitbar не знал, удобная штука))) Да, ребят такой вопрос, как выделить всю оперативную память матлаб? Много где читал что по умолчанию 512 мб стоит, а когда идут тяжелые вычисления он мол всю задействует, но смотря как он считает, я так думаю что 512 и стоит. 0