error message: integral calculation

I want to calculate a formula as followed.

So I wrote a program as below. But in this program, error message" integralCalc/finalInputChecks (line 522). input function must be type of single or double. function hundle found." was displayed.

I tried to use fun1*fun2 in spite of fun, but error messages"function mtimes(input parameter of type 'function_handle' not defined.)" and integralCalc/finalInputChecks (line 314) displayed.

What should I do?

a = 0.2;      %field loss coefficient α
y = 1.3;       % fiber non-linearity coefficient　γ
Ns = 20;        %number of span
Ls = 100;       %span length         
b2 = 20.7;     %dispersion coefficient
%% formula of　ρ

PI = pi;
fun1 = @(z)exp(2*z-2*a*z);
Le = (abs((integral(fun1,0,100))))^2 ;%square of span effective length
fun = @(z,f,f1,f2)exp(1j*4*PI^2*(f1-f)*(f2-f)*b2*z)*exp(2*z-2*a*z);%fun1*fun2
%fun2 = @(f,f1,f2)exp(1j*4*PI^2*(f1-f)*(f2-f)*b2*z); %second element of integration
p = (abs(integral(@(z) fun,0,Ls)))^2/Le%ρ

Thank you.

[Torsten]

So you set g = 1 and beta3 = 0 ?

You didn't give values to f, f1 and f2. If you use "integral", you must give numerical values to all the constants involved except z .

Or do you aim at a symbolic solution ? Then use "int".

Answer 1

[Walter Roberson]

 p = (abs(integral(@(z) fun,0,Ls)))^2/Le%ρ

@(z)fun is an anonymous function that ignores its input argument and returns the anonymous function handle fun. It does not execute fun on the value of z. Just pass fun at that point instead of @(z)fun

[柊介 小山内]

I rewrote program as you said, but another error displayed. 'not enough input parameter' in fun.

error: int>@(z,f,f1,f2)exp(1j*4*PI^2*(f1-f)*(f2-f)*b2*z)*exp(2*z-2*a*z) (行 16)
fun = @(z,f,f1,f2)exp(1j*4*PI^2*(f1-f)*(f2-f)*b2*z)*exp(2*z-2*a*z)%fun1*fun2

error: integralCalc/iterateScalarValued (行 314)
                fx = FUN(t);

error: integralCalc/vadapt (行 132)
            [q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

error: integralCalc (行 75)
        [q,errbnd] = vadapt(@AtoBInvTransform,interval);

error: integral (行 87)
Q = integralCalc(fun,a,b,opstruct);

error: int (行 17)
p = (abs(integral(fun,0,Ls)))^2/Le%ρ

If I write p as follow, error message'parameter or function z don't recognized' displayed.

p = (abs(integral(fun,z,0,Ls)))^2/Le%ρ

What should I do? I'm sorry for spending your time for me.

[Torsten]

You still didn't tell us what f, f1 and f2 are and if you set g = 1 and beta3 = 0 in your original formula.

Knowing this, we can answer your question about the error message.

[柊介 小山内]

f1,f2,f are symbolic parameter and each parameter indicate frequency [THz]. I want to calclate ρ as a function of ,f1,f2,f. I'll integrate ρ respect to f1, f2 after calclating ρ. Finally formula respect to f will be calclated.

l set g(ε)=1 and beta3=0.

Thank you.

[Torsten]

You can't use "integral" in combination with symbolic parameters f1, f2 and f.

syms f1 f2 f a z Ls b2
rho = int(exp(1j*4*sym('pi')^2*(f1-f)*(f2-f)*b2*z)*exp(2*z*(1-a)),z,0,Ls)

rho = 

[Walter Roberson]

The fun that you pass to integral() must be a handle to a function that accepts one parameter.

If the "arrayvalued" option is not used, then the function must return single precision or double precision results exactly the same size as its input (which will not be a fixed size), and for each input value, the result must be the same as-if you had invoked the function on that scalar (so be sure to use element-by-element operations rather than matrix operations.)

If the "arrayvalued" option is used, then the function will be passed a scalar value each time, and must return a single precision or double precision array that is the same size each time.

The function must be deterministic -- no random numbers.

All of this does not mean that the calculation cannot involve symbolic variables -- but the calculation must return single or double precision. For example it would be legal for the calculation to do double(vpasolve(Expression,x))

f1,f2,f being symbolic parameters is a problem if the evaluation of the function cannot get rid of f1, f2, f by the time of return.

If you are trying to do a symbolic integration that returns in terms of f, f1, f2, then you cannot use integral and must instead use int

[柊介 小山内]

thank you for teaching. I can calclate ρ, but I can't integrate ρ respect to f1,f2.

alfa = 0.2;      %field loss coefficient α
gamma = 1.3;       % fiber non-linearity coefficient　γ
Ns = 20;        %number of span
Ls = 100;       %span length         
beta2 = 20.7;     %dispersion coefficient
%% formula of　ρ

PI = pi;
syms z f f1 f2
assume(z>0);Le = (abs((int(exp(-4*alfa*z^2),z,0,Ls))))^2;
assume(f,'real');   
assume(f1,'real');  
assume(f2,'real'); 
Le = (abs((int(exp(-4*alfa*z^2),z,0,Ls))))^2;
p = (abs(int(exp(1j*4*pi^2*(f1-f)*(f2-f)*beta2*z)*exp(2*z*(1-alfa)),z,0,Ls)))^2/Le
%% formula of GNLI 
pint1 = int(p,f1,-0.5,0.5);%integration about f1
pint2 = int(pint1,f2,-0.5,0.5);%integration about f2

error message didn't displayed, but calclation didn't finish.

calclation still continue time goes by 5 minutes.

Is this a problem of my PC(windows10) spec? or something wrong with my code?

notice:

ρ=

(16*abs(a*i/(b*f*f1 + b*f*f2 - b*f1*f2 - b*f^2 + ci) - (exp((f^2*d*i)/g)*exp(-(f*f1*d*i)/g)*exp(-(f*f2*d*i)/g)*exp((f1*f2*d*i)/g)*exp(160)*h*i)/(a*f*f1 + a*f*f2 -a*f1*f2 - a*f^2 + k*i))^2)/(5*pi*erf(40*5^(1/2))^2)

alphabet a～h and k indicate very large constant integer.

I want to calclate formula as below,so I have to integrate ρ respect to f1,f2.( assume all Gwdm=1)

thank you.

[Torsten]

As you can see, "int" cannot determine pint1 and pint2 as analytical expressions.

So I suggest calculating rho and using "integral2" to compute the remaining double integral with respect to f1 and f2.

alfa = 0.2;      %field loss coefficient α
gamma = 1.3;       % fiber non-linearity coefficient　γ
Ns = 20;        %number of span
Ls = 100;       %span length         
beta2 = 20.7;     %dispersion coefficient
%% formula of　ρ

PI = pi;
syms z 
syms f f1 f2 real

assume(z>0);
Le = int(exp(-4*alfa*z^2),z,0,Ls)

Le = 

0.9908

p = int(exp(1j*4*pi^2*(f1-f)*(f2-f)*beta2*z)*exp(2*z*(1-alfa)),z,0,Ls)

p = 

p = real(p)^2+imag(p)^2

p = 

p = p/Le^2;
%% formula of GNLI 
pint1 = int(p,f1,-0.5,0.5)%integration about f1

pint1 = 

pint2 = int(pint1,f2,-0.5,0.5)%integration about f2

pint2 = 

[柊介 小山内]

Thank you. but I want to know pint2 as a function of f. Do you know how to erase f1 and f2? I think integral2 can't use because I want to integrate only f1 and f2, not including f.

[Torsten]

Then - by using integral2 - you get the function G_NLI which is a function of f that can be evaluated by giving f a numerical value.

[柊介 小山内]

if using integral2, How can I specify parameter I want to integrate respect to?

I think integral2 can use only the folmula has two symboric parameter. ρ has three symboric parameter, so I think integral2 can't use.

[Torsten]

fun = @(f1,f2,f) f1.^2+f2.^2+f.^2;
f1min = -1;
f1max = 2;
f2min = 3;
f2max = 6;
G_NLI = @(f)integral2(@(f1,f2)fun(f1,f2,f),f1min,f1max,f2min,f2max);
G_NLI(3)

ans = 279.0000

You will first have to convert rho into a numerical function using "matlabFunction".

[Walter Roberson]

You are attempting to use integral2() to integrate an expression with more than two free variables. It cannot be used for that purpose. Only int() or vpaintergral() can be used in cases where there are more free variables than variables being integrated over (and in such a situation vpaintegral() would just return back the formula unchanged)

What Torsten is showing you is how to create functions that given specific values for the other variables can do numeric integration over particular variables. The result is not a formula: it is a process to do numeric integration.

See also integral3 for numeric integration of order 3, and https://www.mathworks.com/matlabcentral/fileexchange/47919-integraln-m for orders 4, 5, and 6

[柊介 小山内]

Roberson, I used int() at first, but if I use int(),I can't remove f1 and f2 from result NLI. I want to get NLI as a function of f(I have to plot NLI as a function of f ).

PI = pi;
syms z f f1 f2
assume(z>0);Le = (abs((int(exp(2*alfa-2*alfa*z),z,0,Ls))))^2;
assume(f,'real');   
assume(f1,'real');  
assume(f2,'real'); 
Le = ((int(exp(2*alfa-2*alfa*z),z,0,Ls)))^2;
p = (int(exp(1j*4*pi^2*(f1-f)*(f2-f)*beta2*z)*exp(2*z*(1-alfa)),z,0,Ls))^2;
p = real(p)^2*imag(p)^2;
p = p/Le;
%% formula of GNLI
%%p = matlabFunction(p);

pint1 = int(p,f1,-0.5,0.5)%integration about f1
pint2 = int(pint1,f2,-0.5,0.5)%integration about f2
NLI = (16/27)*gamma^2*Le*pint2 % NLI of 1span

if only int and vpaintegral can use, is it impossible to remove f1 and f2?

I'm not good at english, so I 'm sorry if you explained about my question before.

Thank you.

[Torsten]

I want to get NLI as a function of f(I have to plot NLI as a function of f ).

That's exactly what my suggested method does.

[柊介 小山内]

I used integral2 and error message 'input function must be single or double. sym be found.'displayed. If used ht instead of ht(f1,f2,f), sym changed into function handle.

p = (int(exp(1j*4*pi.^2*(f1-f)*(f2-f)*beta2*z)*exp(2*z*(1-alfa)),z,0,Ls))^2;
p = real(p)^2+imag(p)^2;
p = p/Le;
%% formula of GNLI
ht = matlabFunction(p)
pint = @(f)integral2(@(f1,f2)ht(f1,f2,f),-0.5,0.5,-0.5,0.5)
%pint1 = int(p,f1,-0.5,0.5)%integration about f1
%pint2 = int(pint1,f2,-0.5,0.5)%integration about f2
NLI = (16/27)*gamma^2*Le*pint % NLI of 1span

[Torsten]

Are you sure about exp(2*z*(1-alfa)) in the function p ? Shouldn't it be exp(-2*z*(1-alfa)) ? Or exp(-2*alfa*z) as in your graphics ?

alfa = 0.2;      %field loss coefficient α
gamma = 1.3;       % fiber non-linearity coefficient　γ
Ns = 20;        %number of span
Ls = 100;       %span length         
beta2 = 20.7;     %dispersion coefficient
%% formula of　ρ

PI = pi;
syms z 
syms f f1 f2 real

assume(z>0);
Le = int(exp(-4*alfa*z^2),z,0,Ls);
p = int(exp(1j*4*pi^2*(f1-f)*(f2-f)*beta2*z)*exp(2*z*(1-alfa)),z,0,Ls);
p = real(p)^2+imag(p)^2;
p = p/Le^2;
ht = matlabFunction(p)

ht = function_handle with value:
    @(f,f1,f2)(1.0./erf(sqrt(5.0).*4.0e+1).^2.*((imag((exp(f.^2.*8.172032444101988e+4i).*exp(f.*f1.*-8.172032444101988e+4i).*exp(f.*f2.*-8.172032444101988e+4i).*exp(f1.*f2.*8.172032444101988e+4i).*exp(1.6e+2).*7.0368744177664e+14i)./(f.*f1.*5.750556604705831e+17+f.*f2.*5.750556604705831e+17-f1.*f2.*5.750556604705831e+17-f.^2.*5.750556604705831e+17+1.125899906842624e+15i))-((f.*f1.*5.750556604705831e+17+f.*f2.*5.750556604705831e+17-f1.*f2.*5.750556604705831e+17-f.^2.*5.750556604705831e+17).*7.0368744177664e+14)./((f.*f1.*5.750556604705831e+17+f.*f2.*5.750556604705831e+17-f1.*f2.*5.750556604705831e+17-f.^2.*5.750556604705831e+17).^2+1.267650600228229e+30)).^2+(real((exp(f.^2.*8.172032444101988e+4i).*exp(f.*f1.*-8.172032444101988e+4i).*exp(f.*f2.*-8.172032444101988e+4i).*exp(f1.*f2.*8.172032444101988e+4i).*exp(1.6e+2).*7.0368744177664e+14i)./(f.*f1.*5.750556604705831e+17+f.*f2.*5.750556604705831e+17-f1.*f2.*5.750556604705831e+17-f.^2.*5.750556604705831e+17+1.125899906842624e+15i))-7.922816251426434e+29./((f.*f1.*5.750556604705831e+17+f.*f2.*5.750556604705831e+17-f1.*f2.*5.750556604705831e+17-f.^2.*5.750556604705831e+17).^2+1.267650600228229e+30)).^2).*(1.6e+1./5.0))./pi

pint = @(f)integral2(@(f1,f2)ht(f,f1,f2),-0.5,0.5,-0.5,0.5);
NLI = @(f)(16/27)*gamma^2*double(Le)*pint(f) % NLI of 1span

NLI = function_handle with value:
    @(f)(16/27)*gamma^2*double(Le)*pint(f)

NLI(0.1)

ans = 2.2036e+137

[Walter Roberson]

The below code takes several hours to run, mostly the last step. The function handle it produces starts with

@(f)exp(-3.96e+2./5.0).*integral(@(f2)integral(@(f1)exp(3.96e+2./5.0).*1.804851387845415e-35.*(

with the nested calls to integral(). It contains quite a number of very large coefficients beyond 1e250.

The portion I looked at did not have any calls to real() or imag() -- but I can only see a small portion because it is such a long expression.

If you do take this approach I recommend that you use the 'file' option for matlabFunction, and that you specifically ask for optimization to be false when you use the 'file' option (optimization has been producing incorrect results in recent MATLAB releases; I do not know if it is fixed in R2022b.)

This all gives you a NLI as a symbolic expression in f, and gives you NLIF as function handle to process values of f numerically. It is not a closed-form formula: you are unlikely to get a closed form formula unless you do something like taylor series (keeping in mind that taylr series of an exponential is pretty inaccurate unless you restrict the range of inputs a fair bit.)

alfa = 0.2;      %field loss coefficient α
gamma = 1.3;       % fiber non-linearity coefficient　γ
Ns = 20;        %number of span
Ls = 100;       %span length         
beta2 = 20.7;     %dispersion coefficient

PI = sym(pi);
syms z f f1 f2
assume(z>0);
Le = (abs((int(exp(2*alfa-2*alfa*z),z,0,Ls))))^2;
assume(f,'real');   
assume(f1,'real');  
assume(f2,'real'); 
Le = ((int(exp(2*alfa-2*alfa*z),z,0,Ls)))^2;
p = (int(exp(1j*4*pi^2*(f1-f)*(f2-f)*beta2*z)*exp(2*z*(1-alfa)),z,0,Ls))^2;
p = real(p)^2*imag(p)^2;
p = expand(p)/Le;

%% formula of GNLI
%%p = matlabFunction(p);

pint1 = int(p,f1,-0.5,0.5); %integration about f1
pint2 = int(pint1,f2,-0.5,0.5); %integration about f2
NLI = (16/27)*gamma^2*Le*pint2; % NLI of 1span

NLIF = matlabFunction(NLI);

[柊介 小山内]

Thank you Walter and Torsten. If I face new question, I may ask you again.