I am trying solve an integration problem. The result takes around 4-5 minutes to solve. I want to reduce the runtime to max 10seconds. Please help me with the same.

w_p=2*pi/T_p;
syms w;
H_s=5.4;
WD=42.5;
T_p=10.4;
%Generalised Phillip's Constant
alpha=vpa(5*H_s^2*w_p^4*(1-0.287*log(gamma1))/(16*g^2),3);
S_nn=alpha*g^2*w^(-5)*(exp(-(5/4)*(w/w_p)^(-4)))*(gamma1^(exp(-0.5*((w-w_p)/(sigma*w_p))^2)));
syms n;
ss=(w^2)*WD/g;
f=1+0.666*(ss)+0.445*(ss)^2-0.105*(ss)^3+0.272*(ss)^4;
lambda=(2*pi/w)*sqrt(WD*g*(f/(1+ss*f)));
k=2*pi/lambda;
G(w)=w/(sinh(k*WD));
S_UU=G(w)^2*S_nn;
M=vpa(int(S_UU*w^n,0,12.78));
tt=double(M(0));
set(handles.text119,'String',tt);

[Torsten]

And you get a solution without specifying n ?

Best wishes

Torsten.

[Walter Roberson]

Perhaps the integral is with respect to n? However w is also declared symbolic so perhaps the integral is with respect to it. The conciseness of omitting the variable of integration is not worth it when it comes time to debug the problem.

Definitions are missing for g, gamma1, and sigma.

You cannot do effective numeric integration (as implied by the vpa()) when you have free variables; either w or n is a free variable in the integration. However, if the integration happens to be with respect to n, the symbolic integration turns out to be practical.

[Walter Roberson]

Note that your M is not a function, so M(0) would be interpreted as trying to index M at location 0. If you want to evaluate M at n = 0 then use subs(M,n,0)

[Akash_Sharma3008]

Thanks for your responses. I did made a few corrections in the code. I am again posting it here for your easiness. Please go through it. I am facing major problem w.r.t. the time, it is consuming. Accuracy level is all good from my side. Please suggest me some good suggestion that can resolve my query. And I am very much thankful to Mr. Walter Roberson. Sir I am waiting for your response.

T_p=10;
H_s=5;
gamma=0.07;
g=9.81;
WD=40;
w_p=2*pi/T_p;
syms w;
%syms S_nn(w);
sigma=0.07;
format short;
%Generalised Phillip's Constant
alpha=5*H_s^2*w_p^4*(1-0.287*log(gamma))/(16*g^2);
disp(alpha);              %Display- just for the sake of checking
%Spectral Density function 
S_nn(w)=vpa(alpha*g^2*w^(-5)*(exp(-(5/4)*(w/w_p)^(-4)))*(gamma^(exp(-0.5*((w-w_p)/((sigma)*w_p))^2))),4);
disp(S_nn(10));           %Display- just for the sake of checking
syms k;
%syms kk(w);
ss(w)=w^2*WD/g;
f(w)=1+0.666*ss+0.445*ss^2-0.105*ss^3+0.272*ss^4;
lambda(w)=vpa(2*pi/w*sqrt(WD*g*f/(1+ss*f)),4);
disp(lambda(10));        %Display- just for the sake of checking
kk(w)=2*pi/lambda;
disp(kk(1));             %Display- just for the sake of checking 
G(w)=vpa(w/(sinh(kk*WD)),4);
disp(G(10));             %Display- just for the sake of checking
S_UU(w)=vpa((G(w))^2*S_nn(w),4);
disp(S_UU(5));           %Display- just for the sake of checking
M_n=int(w*S_UU(w),0,1);
disp(double(M_n));       %Display- just for the sake of checking

[Walter Roberson]

You define ss as a function in w and then you have

f(w)=1+0.666*ss+0.445*ss^2-0.105*ss^3+0.272*ss^4;

You should be calling ss, not just mentioning its name:

f(w)=1+0.666*ss(w)+0.445*ss(w)^2-0.105*ss(w)^3+0.272*ss(w)^4;

You have similar difficulties in some other locations.

Other than that: as I suggested before, lower your Digits setting. The default is 32 I read recently.

[Akash_Sharma3008]

Thanks for your reply. I have tried both of the method you had suggested. Even after using digits the integration time is 4-5 mins. The same code on Mathcad runs in seconds. Is there any other method that you can suggest? I would be highly grateful to you.

[Walter Roberson]

You could try,

f = w*S_UU(w);
fh = matlabFunction(f,w);
M_n = integrate(fh, 0, 1);

Answer 1

[Walter Roberson]

Provided that the integral is with respect to n, then it is

-(729/58492928) * Pi^4 * exp(-(3125/1827904) * Pi^4/w^4) * gamma1^(exp(-(1/50)  * (-26*w + 5*Pi)^2 /(sigma^2 * Pi^2))) * (-1000 + 287 * ln(gamma1)) * (w^(639/50)-1) / (ln(w) * w^3 * (cosh( (1/2) * w * 85^(1/2) / ((170382840*g^2*w^8 - 1547595*g^3*w^6 + 154326*g^4*w^4 + 5440*g^5*w^2 + 192*g^6) / (14482541400*w^10 - 131545575*g*w^8 + 13117710*g^2*w^6 + 462400*g^3*w^4 + 16320*g^4*w^2 + 384*g^5))^(1/2))^2 - 1))

Assuming that the coefficients should be real and that the results should be real, then the equations have singularities for all non-positive g but not for positive g.

[Akash_Sharma3008]

Thanks Sir for the prompt reply. Sorry I forgot to mention the variables.

g=9.81
sigma=0.07

and gamma1 is obtained from the below mentioned code where phy is the user input. The variable of integration is w and n is just a variable parameter. Kindly help. Thanks

if(phy<=3.6)  
    gamma1=5;
else
    if(3.6<phy<5)
        gamma1=exp((5.75-1.15*phy));
    else
        gamma1=1;
    end
end

[Walter Roberson]

if(3.6<phy<5) will always be true. You want

if (3.6 < phy & phy < 5)

[Walter Roberson]

If you are integrating with respect to w and n does not have a value then you cannot use numeric integration, and there does not happen to be a closed form symbolic integration for that formula.

If your M(0) is intended to mean to evaluate the integral at n = 0, then you should make that substitution before the int() call; the result would be an expression dependent only on w and so it can be integrated numerically.

You might want to decrease DIGITS for faster numeric integration.

I do not think that it is a good idea to use vpa() on the calculation of alpha; the effect on accuracy could be quite large.