Speed up numeric integral

6 visualizzazioni (ultimi 30 giorni)
Guan Hao
Guan Hao il 24 Lug 2024
Modificato: Guan Hao il 25 Lug 2024
Hi,everyone.I got a problem of doing numeric integral on the matrix H (which is at the bottom part of my code).
The code works fine but I would like to speed it up? Am I doing something stupid to make it slow?Thanks for your help.
clear
L=0.1;
section=50;
a=L./2;
b=L./section;
v=3e8;
f1=1.2e9;
f2=4.8e9;
f3=7.5e9;
f4=10e9;
w1=(2.*pi.*f1);
w2=(2.*pi.*f2);
w3=(2.*pi.*f3);
w4=(2.*pi.*f4);
Z01=50;
Z02=75;
a0=(log(Z02./Z01))./(2..*L);
T=a;
Ub=a;
Lb=-a;
ub=a;
lb=-a;
k=5;
syms x y w m
%%% Revised Algorithm
l=0:(2.*L)./(2.*k):L % Approximate sections dvided
% Approximate terms
sec=numel(l)-1;
Z1=[50 55 60 65 70 80 75]; % Set the impedance of the end of each sections
ac=10; % Choose the accuracy of dividing a section
for s=1:1:sec
c=(l(s+1)-l(s))./ac;
app_cos=0.5./ac.*log(Z1(s+1)./Z1(s)).*sin(ac./2.*(w./v).*c).*cos(((ac+1)./2).*(w./v).*c)./sin(((c.*w)./(2.*v)));
app_sin=0.5./ac.*log(Z1(s+1)./Z1(s)).*sin(ac./2.*(w./v).*c).*sin(((ac+1)./2).*(w./v).*c)./sin(((c.*w)./(2.*v)));
amp=0.75;
Ga(s)=(app_cos-(1i.*amp.*app_sin));
end
Ta=(1-(abs(Ga(1:1:end)).^2));
%%% Turn closed loop terms into matrix presentation
for t=1:1:sec
for r=1:1:sec
E(t,r)=Ga(t).*Ga(r);
end
end
E=triu(E,1);
for t=1:1:sec
for r=1:1:sec
D(t,r)=abs(E(r,t));
end
end
C1=sum(D);
% Construct approximate loop terms
for t=1:1:sec
D1(t)=prod(Ta(1:t))./(1+(sum(C1(1:t))));
end
D=[1 D1(2:end)];
Dm=transpose(D).*D;
% Cmn
m=transpose(1:5);
n=1:5;
g=(v.^2.*w.^2.*cos((2.*w.*(- y + x))./v).*cos(20.*sym(pi).*m.*x).*cos(20.*sym(pi).*n.*y) + 10.*sym(pi).*m.*v.^3.*w.*sin(20.*sym(pi).*m.*x).*sin((2.*w.*(- y + x))./v).*cos(20.*sym(pi).*n.*y) - 10.*sym(pi).*n.*v.^3.*w.*sin(20.*sym(pi).*n.*y).*sin((2.*w.*(- y + x))./v).*cos(20.*sym(pi).*m.*x) + 100.*m.*n.*v.^4.*sym(pi).^2.*sin(20.*sym(pi).*m.*x).*sin(20.*sym(pi).*n.*y).*cos((2.*w.*(- y + x))./v))./(4.*(w.^2 - 100.*sym(pi).^2.*m.^2.*v.^2).*(w.^2 - 100.*sym(pi).^2.*n.^2.*v.^2));
for p=1:numel(l)-1
x1=l(p);
x2=l(p+1);
for q=1:numel(l)-1
y1=l(q);
y2=l(p+1);
Ht{p,q}=matlabFunction(((subs(g,[x y],[x2 y2])-subs(g,[x y],[x1 y1]))*Dm(p,q)));
end
end
Hi=Ht{1};
for m=2:1:numel(Ht)
Hi=Hi+vpa(Ht{m},3);
end
H=integral(Hi,w1,w2,"ArrayValued",true)
  4 Commenti
Torsten
Torsten il 24 Lug 2024
Modificato: Torsten il 24 Lug 2024
Do you see your mistake ?
syms x y
x1 = 1;
x2 = 2;
y1 = 4;
y2 = 6;
g = 4*x*y;
G = int(int(g,x),y)
G = 
int(int(g,x,x1,x2),y,y1,y2)
ans = 
60
int(int(g,y,y1,y2),x,x1,x2)
ans = 
60
subs(G,[x y],[x2,y2])-subs(G,[x y],[x1,y1])
ans = 
128
Guan Hao
Guan Hao il 24 Lug 2024
Modificato: Guan Hao il 25 Lug 2024
@Torsten Oh.Thanks for pointing that out.

Accedi per commentare.

Risposte (0)

Categorie

Scopri di più su MATLAB in Help Center e File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by