Integration for a range of values with variable limits

Hi all,

I need to solve the following equation for each individual value of f where f=0:0.01:10:

The equation is for MSAR in the code below, with the intention of getting a range of RMSAR values for a given f.

K1=274350; C1=11040; K2=414000; C2=16000;
M=13200; L1=2.34; L2=2.885; K=2.310;
I=M*(K^2);
a1=46.85*(10^-4); b1=0.19;

syms f

w=f*2*pi;
B11=K1+K2-M*(w^2)+1i*w*(C1+C2);
B12=K2*L2-K1*L1+1i*w*(-C1*L1+C2*L2);
B13=-K1-1i*w*C1;
B21=K2*L2-K1*L1+1i*w*(-C1*L1+C2*L2);
B22=K1*(L1^2)+K2*(L2^2)-I*(w^2)+1i*w*(C1*L1^2+C2*L2^2);
B23=K1*L1+1i*w*C1*L1;
B31=0; B32=0; B33=1;
S1=K2+1i*w*C2;
S2=K2*L2+1i*w*C2*L2;
S3=0.975;

B(1,1)=B11; B(1,2)=B12; B(1,3)=B13;
B(2,1)=B21; B(2,2)=B22; B(2,3)=B23;
B(3,1)=B31; B(3,2)=B32; B(3,3)=B33;
S(1,1)=S1; S(2,1)=S2; S(3,1)=S3;

T=B\S;
hif=T(1,:)';
G1=a1*exp((-b1)*f);

MSARn=@(f)((abs(hif))^2).^(f^4)*G1;

iMSARn=int(MSARn,0.89*f,1.12*f);

MSAR=((2*pi)^4)*iMSARn;

RMSAR=MSAR.^0.5;

When I try running the code, it's incredibly slow. I left it going over night and it still wasn't done. So I need a more practical method as this is only a section of a larger code. Am I correct in thinking that only 'int' accepts variable integral limits? I've tried 'quad' in a loop but it still sees the limits as non-scalar:

for k=1:1:1001;
    f(k)=(k-1)/100;
    iMSARn(k)=quad(@PSD_Opt_int_loop,[0.89*f(k),1.21*f(k)]);
end

Thanks in advance.

Answer 1

[Walter Roberson]

      iMSARn(k)=quad(@PSD_Opt_int_loop, 0.89*f(k), 1.21*f(k));

Notice no []

[Keelan Toal]

That solves that problem, thanks.

However, I'm having further issues. The code runs entirely but my results are way out. Can you see any obvious issues with my code? Assuming my constants/equations are correct.

function MSARn = Func_PSD2 (f)
K1=274350; C1=11040; K2=414000; C2=16000;
M=13200; L1=2.34; L2=2.885; K=2.310;
I=M*(K^2);
a1=46.85*(10^-4); b1=0.19;

w=f*2*pi;
B11=K1+K2-M*(w.^2)+1i*w*(C1+C2);
B12=K2*L2-K1*L1+1i*w*(-C1*L1+C2*L2);
B13=-K1-1i*w*C1;
B21=K2*L2-K1*L1+1i*w*(-C1*L1+C2*L2);
B22=K1*(L1^2)+K2*(L2^2)-I*(w.^2)+1i*w*(C1*(L1^2)+C2*(L2^2));
B23=K1*L1+1i*w*C1*L1;
B31=0; B32=0; B33=1;
S1=K2+1i*w*C2;
S2=K2*L2+1i*w*C2*L2;
S3=0.975;

b11=size(B11,2);
B11p=permute(reshape(B11',1,b11,[]),[1 3 2]);
B12p=permute(reshape(B12',1,b11,[]),[1 3 2]);
B13p=permute(reshape(B13',1,b11,[]),[1 3 2]);
B21p=permute(reshape(B21',1,b11,[]),[1 3 2]);
B22p=permute(reshape(B22',1,b11,[]),[1 3 2]);
B23p=permute(reshape(B23',1,b11,[]),[1 3 2]);
S1p=permute(reshape(S1',1,b11,[]),[1 3 2]);
S2p=permute(reshape(S2',1,b11,[]),[1 3 2]);

B=zeros(3,3,b11);

B(1,1,:)=B11p; B(1,2,:)=B12p; B(1,3,:)=B13p;
B(2,1,:)=B21p; B(2,2,:)=B22p; B(2,3,:)=B23p;
B(3,1,:)=B31; B(3,2,:)=B32; B(3,3,:)=B33;
S(1,1,:)=S1p; S(2,1,:)=S2p; S(3,1,:)=S3;

T(:,:,1)=abs(B(:,:,1))\abs(S(:,:,1));
T(:,:,2)=abs(B(:,:,2))\abs(S(:,:,2));

if size(B,3)>2
    T(:,:,3)=abs(B(:,:,3))\abs(S(:,:,3));
end
if size(B,3)>3
    T(:,:,4)=abs(B(:,:,4))\abs(S(:,:,4));
end
if size(B,3)>4
    T(:,:,5)=abs(B(:,:,5))\abs(S(:,:,5));
end
if size(B,3)>5
    T(:,:,6)=abs(B(:,:,6))\abs(S(:,:,6));
end
if size(B,3)>6
    T(:,:,7)=abs(B(:,:,7))\abs(S(:,:,7));
end

hif=((squeeze(T(1,:,:)))');
G1=a1*exp((-b1)*f);

MSARn=(hif.^2).*(f.^4).*G1;
end

and

clear
clc
f=zeros(1,1001); iMSARn=zeros(1,1001);
for k=1:1:1001;
    f(k)=(k)/100;
    iMSARn(k)=quad(@Func_PSD2,0.89*f(k),1.12*f(k));
end
MSAR=((2*pi)^4)*iMSARn;
RMSAR=MSAR.^0.5;
plot(f,RMSAR)

I'm aiming for the graph labelled 4:

But I'm getting the following:

Any suggestions?