Sum of difference of integrals

Hello,

I've tried to visualise the value of a function F below (that it's the sum of the difference between two integrals, between 10e3 et 30.1e6 ), but the values of the function F when running the programm is "NaN" . Could you please tell me when I am wrong?

The program is below . Thank you enourmously !

%% Sth
z= 76.05*1e-2;
delta=6*1e-3;
c = 3e8;
muo=4*pi*10^-7;
fr= [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur = 1;
w=2*pi*fr;
lambda = c./fr;
gamao= 1j*w/c;
s=1;
gama= sqrt(1j.*2.*pi.*fr.*muo.*mur.*s);
too= sqrt(lambda.^2 -gamao.^2);
to= sqrt(lambda.^2- gama.^2);
K= ((to./too+mur).^2-(to./too-mur).^2.*exp(-2.*to.*delta)).^-1;
R = 16.05*1e-2;
J = besselj(1,lambda*R);
L = lambda.*lambda.*((too).^-1).*J.*exp(-too.*z);
G= int(L,sym('lambda'),0,inf);
M = K.*(lambda.^2).*(to).*(too.^-2).*J.*exp(-too*z-(to-too)).* delta;
H= int(K.*(lambda.^2).*(to).*(too.^-2).*J.*exp(-too*z-(to-too)),sym('lambda'),0,inf);
Sth= abs((1/4*mur).*(G./H));
% %% Se
Se = [1, 2, 3, 6, 9 ,10 , 12 ,14 ,18, 22 ,24, 28];
%% F calculation
F= symsum(abs(Se-Sth).^2,f,10e3,30.1e6);

Answer 1

[Walter Roberson]

G= int(L,sym('lambda'),0,inf);

the sym('lambda') there has nothing to do with the lambda in the previous lines, which are numeric lambda.

Your L is numeric, and int() of a constant is the constant times the difference in bounds.

You symsum() with respect to f, but the inputs are an array of numeric values, and you have no variable named f to sum over.

symsum() should typically only be used when you need to generate a formula from the summation, because you have an indefinite bound. If you have a vector of values, then you should sum() the vector instead. If you have definite upper and lower bound and a formula that you are symsum() then typically you should instead subs() lower:upper into the formula and sum() that.

[Salwa Ben Mbarek]

Hello Mr. Roberson,

Thank you very much for your help and your explanations.

I was not clear sorry , the G and H integrals are supposed to be depending on lambda (the variable is lambda), here is the formula in the picture :

Is it no more a constant ? when I write the whole expression in G :

G= int(lambda.*lambda.*((too).^-1).*J.*exp(-too.*z),lambda,0,inf);

I get : Check for missing argument or incorrect argument data type in call to function 'int'.

Can you please help me? Thank you !

[Walter Roberson]

you have

c = 3e8;
fr= [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
lambda = c./fr;

Those define lambda as being a numeric vector of 12 values. You try to integrate with respect to lambda but you cannot integrate with respect to a number.

You need to decide whether you are using numeric formulas without integration, or if you are doing useless numeric integration (integration of numeric values instead of a formula), or if you are doing useless symbolic integration, or if you used the wrong formulas and meant to have a symbolic variable in there somewhere.

I would suggest that you make fr symbolic and integrate with respect to fr, and afterwards you can substitute the list of numeric values for fr.

[Salwa Ben Mbarek]

Dear sir,

Thank you so much for your help . I have defined fr symbolic as you suggested so I can substitude the list of numeric values for fr, s symbolic (because I will perform an optimization on the value of s in an other program), but the calculation took hours to be simulated (the matlab status remains busy) . Can you please tell me if there's any update on Matlab or how can I reduce time? Thank you so much.

Here is the program :

%

% Sth
z= 76.05*1e-2;
delta=6*1e-3;
c = 3e8;
muo=4*pi*10^-7;
fr= [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur = 1;
w=2*pi*fr;
lambda = c./fr;
gamao= 1j*w/c;
sym s;
gama= sqrt(1j.*2.*pi.*fr.*muo.*mur.*s);
too= sqrt(lambda.^2 -gamao.^2);
to= sqrt(lambda.^2- gama.^2);
K= ((to./too+mur).^2-(to./too-mur).^2.*exp(-2.*to.*delta)).^-1;
R = 16.05*1e-2;
J = besselj(1,lambda*R);
L = lambda.*lambda.*((too).^-1).*J.*exp(-too.*z);
sym lambda
G= int(L,lambda,0,inf);
M = K.*(lambda.^2).*(to).*(too.^-2).*J.*exp(-too*z-(to-too)).* delta;
H= int(M,0,inf);
o = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
p= c ./ o;
N= subs(G,lambda,p);
Q= subs(H,lambda,p);
SEHth= abs((1/4*mur).*(N./Q));
% %% Se
Se = [1, 2, 3, 6, 9 ,10 , 12 ,14 ,18, 22 ,24, 28];

[Walter Roberson]

I annotate all the variables so you could follow which are scalars or vectors and which are symbolic (and when)

In the below, _Ns indicates Numeric Scalar, _Nv indicates Numeric Vector, _Ss indicates Symbolic Scalar, and _Sv indicates Symbolic Vector.

You need to re-think how you use lambda.

% Sth
z_Nv = 76.05*1e-2;
delta_Ns = 6*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;

%in this section, lambda is a numeric vector with 12 elements
lambda_Nv = c_Ns ./ fr_Nv;

gamao_Nv = 1j.*w_Nv./c_Ns;

%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns .* s_Ss);  %uses s
too_Nv = sqrt(lambda_Nv.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Nv.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Nv + mur_Ns).^2 - (to_Sv./too_Nv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Nv = besselj(1, lambda_Nv .* R_Ns);                                %does not use s
L_Nv = lambda_Nv .* lambda_Nv .* ((too_Nv).^-1) .* J_Nv .* exp(-too_Nv .* z_Nv);   %does not use s
%pay attention to that: L is purely numeric operations with no symbolic s, and lambda is still numeric

%in this section, lambda is a symbolic scalar
lambda_Ss = sym('lambda');

%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = int(L_Nv, lambda_Ss, 0, inf);
%and numeric constant integrated over a range is difference in range times the constant,
%and is 0 where the constant is zero, and otherwise is infinity times sign of the constant
%so G_Sv is full of 0 and infinity

M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Nv.^-2) .* J_Nv .* exp(-too_Nv .* z_Nv - (to_Sv - too_Nv)) .* delta_Ns;
H_Sv = vpaintegral(M_Sv, s, 0, inf);
o_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
p_Nv = c_Ns ./ o_Nv;
N_Sv = subs(G_Sv, lambda_Ss, p_Nv);    %remember we provided G_Sv is all 0 and inf
Q_Sv = subs(H_Sv, lambda_Ss, p_Nv);    %large parts of this one are very small like 1e-9000
SEHth_Sv = abs((1/4 .* mur_Ns) .* (N_Sv ./ Q_Sv));   %N_Sv is all 0 and inf so that controls the outcome
% %% Se
Se_Nv = [1, 2, 3, 6, 9 ,10 , 12 ,14 ,18, 22 ,24, 28];

[Salwa Ben Mbarek]

Dear sir,

Thank you so much for your help! the code is now fast thank you ! I chose to put lambda symbolic in the whole code and substitute at the end. I cannot find the numerical values of N_Sv and Q_Sv (lambda is not replaced by c/.fr) . When running the program, Matlab keeps "lambda" in the expressions of N_Sv and Q_Sv. Here below is the program (I add comments as you did on your previous reply) . Could you please tell me why the value of lambda is not remplaced at the end ? Thank you enormously for all your help !

% Sth
z_Ns = 76.05*1e-2;
delta_Ns = 6*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns .* s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
lambda_Ss = sym('lambda');
G_Sv = int(L_Sv,0, inf);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf);
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
p_Nv = c_Ns ./ fr_Nv;
lambda_Ss = sym('lambda');
N_Sv = subs(G_Sv, lambda_Ss, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
Q_Sv = subs(H_Sv,lambda_Ss , p_Nv);    %Same substitution problem
SEHth_Sv = abs((1/4 .* mur_Ns) .* (N_Sv ./ Q_Sv));   %N_Sv is all 0 and inf so that controls the outcome
% %% Se
Se_Nv = [1, 2, 3, 6, 9 ,10 , 12 ,14 ,18, 22 ,24, 28];

[Walter Roberson]

You should re-check your L_Sv. At the moment the expressions only involve lambda, and when you integrate those with respect to lambda over fixed finite bounds, you are going to get numeric results independent of any variable for G_Sv. But later you expect to substitute lambda values in for G_Sv, which isn't going to work because they are independent of lambda.

% Sth
z_Ns = 76.05*1e-2;
delta_Ns = 6*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns .* s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf);
p_Nv = c_Ns ./ fr_Nv;
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %Same substitution problem
SEHth_Sv = abs((1/4 .* mur_Ns) .* (N_Sv ./ Q_Sv));   %N_Sv is all 0 and inf so that controls the outcome
% %% Se
Se_Nv = [1, 2, 3, 6, 9 ,10 , 12 ,14 ,18, 22 ,24, 28];

output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);

[Salwa Ben Mbarek]

Hello Sir,

Thank you so much for your help, I re-checked the value of L_Sv and you're right it's normal that the values are independent of lambda .

For H_Sv and Q_Sv, l still don't know why Matlab does not calculate the integral, normally we should have also results independent of lambda ,maybe dependent of s or any other variable but not lambda? as you explained we have boundaries . Could you please help me ? Thanks a lot.

Here's the programm :

% Sth
z_Ns = 76.05*1e-2;
delta_Ns = 6*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns .* s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf); % The results depends on lambda 
%%lambda_Ss= c_Ns./fr_Nv;
p_Nv = c_Ns ./ fr_Nv;
%N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, lambda_Ss);    %Cannot substitute and replace the lambda by o and I don't know why
%%Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, lambda_Ss);
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %lambda is in the expression
SEHth_Sv = abs((1/(4 .* mur_Ns)) .* (N_Sv ./ Q_Sv));
% %% Se
Se_Nv = [1, 2, 3, 6, 9 ,10 , 12 ,14 ,18, 22 ,24, 28];
output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);

[Walter Roberson]

Your M_Sv is a matrix of integrals of complex values with a complicated formula. It is not surprising that MATLAB is not able to find a closed-form result.

[Salwa Ben Mbarek]

Dear Sir,

Thank you for your response .

I want to minimise using the interior point method the function Fc below in matlab, to obtain the value of s, knowing that lambda is always c/fr (I've tried to substitute its value in Fc, so the only variable left is s, but I could not). The Fc is calculated through the formula below. I have an error message while running the program "Unrecognized function or variable 'Fc' ". Can you please help me?

Below is the formula and the matlab code. I also tried to use directly the optimization tool at apps but I have the same message error: Optimization running.Error running optimization.Undefined function 'Fc' for input arguments of type 'double'.

Thank you a lot !

Here is the code :

% Sth
z_Ns = 76.05*1e-2;
delta_Ns = 6*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns .* s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf);
%%lambda_Ss= c_Ns./fr_Nv;
p_Nv = c_Ns ./ fr_Nv;
%N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, lambda_Ss);    %Cannot substitute and replace the lambda by o and I don't know why
%%Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, lambda_Ss);
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %Same substitution problem
SEHth_Sv = abs((1/(4 .* mur_Ns)) .* (N_Sv ./ Q_Sv));
% %% Se
Se_Nv = [1, 2, 3, 6, 9 ,10 , 12 ,14 ,18, 22 ,24, 28];
output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);
X= abs((Se_Nv - SEHth_Sv).^2);
F=sum(subs(X,'lambda_Ss',p_Nv));
%% objective harmonization
%Using fmincon command to solve The Robust Reliability Design
obj=@F;
A=[]; 
b=[];
s_Ss0=10e-22;
lb=10e-22;
ub=10e8;
Aeq = [];
beq = [];
nonlcon=@(s_Ss)constrain(s_Ss,fr_Nv,muo_Ns,mur_Ns,gama_Sv);
fmin = fmincon(obj,s_Ss0,A,b,Aeq,beq,lb,ub,nonlcon);

the constraint function is :

% create file constrain.m for nonlinear constraints
function [c,ceq] = constrain(s_Ss,fr_Nv,muo_Ns,mur_Ns,gama_Sv)
  c = sqrt(1j.*2.*pi.*fr_Nv.*muo_Ns.*mur_Ns.*s_Ss);
  ceq= sqrt(lambda_Ss.^2- gama_Sv.^2);

[Walter Roberson]

Your code uses obj=@F where F is a symbolic expression, and attempts to fmincon that function. You cannot minimize a symbolic expression using fmincon. See matlabFunction and be sure to use the vars option passing in a cell array containing a vector that contains the variables to be optimized over.

However your F is constructed using abs(), and abs is not differentiateable at 0. fmincon can only process differentiateable functions. Rewrite your abs()^2 to real() squared plus imag() squared

[Salwa Ben Mbarek]

Dear Sir,

Thank you ! I've tried to do what you suggested. I have an error message for invalid text character. Can you please tell me where I am wrong? Thank you so much. The error messages and the codes are below:

"Warning: Function '_range' not verified to be a valid MATLAB function.

Error using symengine

Error: Invalid text character. Check for unsupported symbol, invisible character, or pasting of non-ASCII characters.

Error in symengine

Error in sym/matlabFunction (line 190)

g = symengine('makeFhandle',varnames,body);

Error in Substi_inch_25_05 (line 45)

Fc= matlabFunction(F,'Vars',{s_Ss});"

Here is the code :

% Sth
z_Ns = 76.05*1e-2;
delta_Ns = 6*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns .* s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf);
%%lambda_Ss= c_Ns./fr_Nv;
p_Nv = c_Ns ./ fr_Nv;
%N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, lambda_Ss);    %Cannot substitute and replace the lambda by o and I don't know why
%%Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, lambda_Ss);
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %Same substitution problem
SEHth_Sv = abs((1/(4 .* mur_Ns)) .* (N_Sv ./ Q_Sv));
% %% Se
Se_Nv = [1, 2, 3, 6, 9 ,10 , 12 ,14 ,18, 22 ,24, 28]
output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);
X= abs((Se_Nv - SEHth_Sv).^2);
X= (real(Se_Nv - SEHth_Sv)).^2+(imag(Se_Nv - SEHth_Sv)).^2;
F=sum(subs(X,'lambda_Ss',p_Nv));
%% objective harmonization
Fc= matlabFunction(F,'Vars',{s_Ss});
obj=@Fc;
A=[]; 
b=[];
s_Ss0=10e-22;
lb=10e-22;
ub=10e8;
Aeq = [];
beq = [];
nonlcon=@(s_Ss)constrain(s_Ss,fr_Nv,muo_Ns,mur_Ns,gama_Sv);
fmin = fmincon(obj,s_Ss0,A,b,Aeq,beq,lb,ub,nonlcon);

[Walter Roberson]

matlabFunction cannot process vpaintegral()

[Walter Roberson]

The below will get you much closer.

Some useless steps have been skipped, and some performance problems have been worked around, and some MATLAB bugs have been worked around.

However... You need to fix your nonlinear constraint function.

• Your nonlinear constraint function tries to use lambda, but lambda has no role at that point. The objective function is strictly based on s
• Your c computes values that are certain to be complex-valued given your upper and lower bounds on s (in particular, the fact that is is real-valued at all). The automatic comparison of those complex values to 0 is probably not going to do what you want
• Your ceq computes values that involve the square of complex quantities that you can be certain are not the square roots of a negative real value -- that is, the square of the quantities is going to give you a complex quantity. Those are then subtracted from the undefined lambda, squared. It is unlikely that the result will exactly balance to be a positive real number suitable for comparison to 0 for equality constraints. The algorithm expects a real number to be returned -- the algorithm uses the sign() of the returned value to try to figure out how to improve the input parameter so that the equality is met. If it does not just fall over, it is unlikely to do what you want.

Also, my tests appear to show that for s values up to about 1e-10, that numerically the result of F is exactly the same, about 281392.545038894 and that the value rises only slowly until about 1e-5. By about 0.25 it has started to take off rapidly.

The point being that the minimum occurs at the lowest permitted (positive) value that satisfies the nonlinear constraints.

... And, knowing that, you can skip a lot of calculations and only do as much as you need to search for locations that satisfy the nonlinear constraints.

% Sth
z_Ns = 76.05*1e-2;
delta_Ns = 6*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns) .* sqrt(s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf);
G_Nv = double(G_Sv);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;

%the integral for H_Sv is going to make it into the final expression. However, matlabFunction cannot work
%with vpaintegral() without falling over, and int() takes a long time to figure out that it is not going to
%be able to find an answer. In the interests of performance, code vpaintegral() here, and we will use
%a hack to swap it out to int() later.

H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf);   %this takes a long time to figure out it cannot do anything.

%%lambda_Ss= c_Ns./fr_Nv;
p_Nv = c_Ns ./ fr_Nv;
%N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, lambda_Ss);    %Cannot substitute and replace the lambda by o and I don't know why
%%Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, lambda_Ss);
%{
G_Sv / G_Nv are independent of any variables so no point doing the subs()
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Nv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
%}
N_Nv = G_Nv;

%{
H_Sv is independent of lambda, so there is no point doing the subs()
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %Same substitution problem
%}
Q_Sv = H_Sv;

temp_Sv = (1/(4 .* mur_Ns)) .* (N_Nv ./ Q_Sv);

SEHth_Sv = sqrt(real(temp_Sv).^2 + imag(temp_Sv).^2);   %avoid using abs, it is not continuously differentiable
% %% Se
Se_Nv = [1, 2, 3, 6, 9 ,10 , 12 ,14 ,18, 22 ,24, 28];

%{
output_Sv is not used so no point computing it
output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);
%}
%X= abs((Se_Nv - SEHth_Sv).^2);     %abs is not differentiable, use the smooth function instead
%{
SEHth_Sv is already abs() so we do not have to worry about imaginary
%X = (real(Se_Nv - SEHth_Sv)).^2+(imag(Se_Nv - SEHth_Sv)).^2;
%}
X = (Se_Nv - SEHth_Sv).^2;

%{
X is independent of lambda, so no point doing the subs
F=sum(subs(X,lambda_Ss,p_Nv));
%}
F = sum(X);
Fint = mapSymType(F, 'vpaintegral', @(V) feval(symengine, 'hold(int)', feval(symengine, 'op', V, 1), 'lambda = 0..infinity'));

%% objective harmonization
Fc = matlabFunction(Fint,'Vars',{s_Ss});   %if you send F to file then optimize must be off
gama_fun = matlabFunction(gama_Sv, 'Vars', {s_Ss}, 'File', 'gama_fun.m', 'optimize', true);
obj = Fc;
A=[]; 
b=[];
s_Ss0=10e-22;
lb=10e-22;
ub=10e8;
Aeq = [];
beq = [];
nonlcon = @(s_Ss) constrain(s_Ss, fr_Nv, muo_Ns, mur_Ns, gama_fun);
fmin = fmincon(obj, s_Ss0, A, b, Aeq, beq, lb, ub, nonlcon);

function [c,ceq] = constrain(s_Ss, fr_Nv, muo_Ns, mur_Ns, gama_fun)
  c = sqrt(1j.*2.*pi.*fr_Nv.*muo_Ns.*mur_Ns.*s_Ss);
  ceq = sqrt(lambda_Ss.^2 - gama_fun(s_Ss).^2);

[Salwa Ben Mbarek]

Dear sir,

First, thank you enormously for your time, your help and your kidness . You helped me a lot !! Many thanks to you !

Thank you also for your explanations, I'm learning so much from you .

Concerning the constraint function, in fact I don't really need a one, I just want to know the value of s_Ss in the frequency range [10kHz- 30 MHz]. The only condition that I have is that gama_Sv should contain s_Ss in : gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns) .* sqrt(s_Ss) , but that is already mentionned in the main program.

Sir, when we sum X , is it always with respect to the frequency range (fmin and fmax) in the program that you corrected ?

Here is the program, I put [] in the constraints and displayed the final result in the program below that gaves me "11039591.337" because I changed the value of delta_Ns . If you find anything incorrect please correct me .Thank you very much !

% Sth
z_Ns = 76.05*1e-2;
delta_Ns = 1*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns) .* sqrt(s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf);
G_Nv = double(G_Sv);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
%the integral for H_Sv is going to make it into the final expression. However, matlabFunction cannot work
%with vpaintegral() without falling over, and int() takes a long time to figure out that it is not going to
%be able to find an answer. In the interests of performance, code vpaintegral() here, and we will use
%a hack to swap it out to int() later.
H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf);   %this takes a long time to figure out it cannot do anything.
%%lambda_Ss= c_Ns./fr_Nv;
p_Nv = c_Ns ./ fr_Nv;
%N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, lambda_Ss);    %Cannot substitute and replace the lambda by o and I don't know why
%%Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, lambda_Ss);
%{
G_Sv / G_Nv are independent of any variables so no point doing the subs()
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Nv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
%}
N_Nv = G_Nv;
%{
H_Sv is independent of lambda, so there is no point doing the subs()
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %Same substitution problem
%}
Q_Sv = H_Sv;
temp_Sv = (1/(4 .* mur_Ns)) .* (N_Nv ./ Q_Sv);
SEHth_Sv = sqrt(real(temp_Sv).^2 + imag(temp_Sv).^2);   %avoid using abs, it is not continuously differentiable
% %% Se
Se_Nv = [12, 19, 21, 26, 29 ,34 , 38 ,39 ,41, 42 ,41, 42];
%{
output_Sv is not used so no point computing it
output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);
%}
%X= abs((Se_Nv - SEHth_Sv).^2);     %abs is not differentiable, use the smooth function instead
%{
SEHth_Sv is already abs() so we do not have to worry about imaginary
%X = (real(Se_Nv - SEHth_Sv)).^2+(imag(Se_Nv - SEHth_Sv)).^2;
%}
X = (Se_Nv - SEHth_Sv).^2;
%{
X is independent of lambda, so no point doing the subs
F=sum(subs(X,lambda_Ss,p_Nv));
%}
F = sum(X);
Fint = mapSymType(F, 'vpaintegral', @(V) feval(symengine, 'hold(int)', feval(symengine, 'op', V, 1), 'lambda = 0..infinity'));
%% objective harmonization
Fc = matlabFunction(Fint,'Vars',{s_Ss});   %if you send F to file then optimize must be off
gama_fun = matlabFunction(gama_Sv, 'Vars', {s_Ss}, 'File', 'gama_fun.m', 'optimize', true);
obj = Fc;
A=[]; 
b=[];
s_Ss0=10e-22;
lb=10e-22;
ub=10e8;
Aeq = [];
beq = [];
% %nonlcon = @(s_Ss) constrain(s_Ss, fr_Nv, muo_Ns, mur_Ns, gama_fun);
nonlcon= @(s_Ss) constrain(s_Ss);
% optimize with fmincon
%[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN]
%%fmin= fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS);
%%fmin = fmincon(obj, s_Ss0, A, b, Aeq, beq, lb, ub, nonlcon);
s_Ss = fmincon(obj,s_Ss0,A,b,Aeq,beq,lb,ub,nonlcon);

% show final objective
disp(['Final Objective: ' num2str(obj(s_Ss))])

constrain function : (file constrain.m)
function [c,ceq] = constrain(s_Ss)
c = [];
ceq = [];
end

[Salwa Ben Mbarek]

Dear sir,

I have changed the value of lambda when we sum X because I want have lambda (c/f) in the frequency range fr.

Here is the new program:

% Sth
z_Ns = 92.18*1e-2;
delta_Ns = 0.8*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
e= c_Ns ./ fr_Nv;
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns) .* sqrt(s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf);
G_Nv = double(G_Sv);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
%the integral for H_Sv is going to make it into the final expression. However, matlabFunction cannot work
%with vpaintegral() without falling over, and int() takes a long time to figure out that it is not going to
%be able to find an answer. In the interests of performance, code vpaintegral() here, and we will use
%a hack to swap it out to int() later.
H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf);   %this takes a long time to figure out it cannot do anything.
%%lambda_Ss= c_Ns./fr_Nv;
p_Nv = c_Ns ./ fr_Nv;
%N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, lambda_Ss);    %Cannot substitute and replace the lambda by o and I don't know why
%%Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, lambda_Ss);
%{
G_Sv / G_Nv are independent of any variables so no point doing the subs()
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Nv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
%}
N_Nv = G_Nv;
%{
H_Sv is independent of lambda, so there is no point doing the subs()
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %Same substitution problem
%}
Q_Sv = H_Sv;
temp_Sv = (1/(4 .* mur_Ns)) .* (N_Nv ./ Q_Sv);
SEHth_Sv = sqrt(real(temp_Sv).^2 + imag(temp_Sv).^2);   %avoid using abs, it is not continuously differentiable
% %% Se
Se_Nv = [12, 19, 21, 26, 29 ,34 , 38 ,39 ,41, 42 ,41, 42];
%{
output_Sv is not used so no point computing it
output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);
%}
%X= abs((Se_Nv - SEHth_Sv).^2);     %abs is not differentiable, use the smooth function instead
%{
SEHth_Sv is already abs() so we do not have to worry about imaginary
%X = (real(Se_Nv - SEHth_Sv)).^2+(imag(Se_Nv - SEHth_Sv)).^2;
%}
X = (Se_Nv - SEHth_Sv).^2;
%{
X is independent of lambda, so no point doing the subs
F=sum(subs(X,lambda_Ss,p_Nv));
%}
F = sum(X);
Fint = mapSymType(F, 'vpaintegral', @(V) feval(symengine, 'hold(int)', feval(symengine, 'op', V, 1), 'lambda= 3.0000..0.0010'));
%% objective harmonization
Fc = matlabFunction(Fint,'Vars',{s_Ss});   %if you send F to file then optimize must be off
gama_fun = matlabFunction(gama_Sv, 'Vars', {s_Ss}, 'File', 'gama_fun.m', 'optimize', true);
obj = Fc;
A=[]; 
b=[];
s_Ss0=10e-22;
lb=10e-22;
ub=10e8;
Aeq = [];
beq = [];
% %nonlcon = @(s_Ss) constrain(s_Ss, fr_Nv, muo_Ns, mur_Ns, gama_fun);
nonlcon= @(s_Ss) constrain(s_Ss);
% optimize with fmincon
%[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN]
%%fmin= fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS);
%%fmin = fmincon(obj, s_Ss0, A, b, Aeq, beq, lb, ub, nonlcon);
s_Ss = fmincon(obj,s_Ss0,A,b,Aeq,beq,lb,ub,nonlcon);
% show final objective
disp(['Final Objective: ' num2str(obj(s_Ss))])
constrain function : (file constrain.m)
function [c,ceq] = constrain(s_Ss)
c = [];
ceq = [];
end

I found a value of s of : 61960391.5578, do you think that the program is right ?

Thank you a lot !

[Walter Roberson]

No, you are running into precision problems. If you try

obj(1e-21)

you will get a value that is slightly below your calculated 61960391.5578: in particular 61960391.5577984 instead of 61960391.5577992 . But you get that 61960391.5577984 value for all values below about 1e-10 .

Smaller inputs lead to smaller outputs in theory, but the difference is smaller than you get measure using double precision.

61960391.557797163359587255798857649933351008580208 : Fint at the location you find using fmincon
61960391.557796377343160728478818104132535439753431 : Fint at 1e-18
61960391.557796377343160728478817973411834316859228 : Fint at 1e-21 (your lower bound)
61960391.557796377343160728478817973411703596027384 : Fint at 0

[Salwa Ben Mbarek]

Thank you so much for your answer. Does it mean that I should keep the initial algorithm ? when :

Fint = mapSymType(F, 'vpaintegral', @(V) feval(symengine, 'hold(int)', feval(symengine, 'op', V, 1), 'lambda = 0..infinity'));  

the final objective is 17427189.0848. Thank you so much !

Here are some tests :

obj(1e-22) =1.7427e+07
obj(1e8) =Inf
obj(0)=1.7427e+07
obj(1e3) =6.9310e+297

The program is :

% Sth
z_Ns = 92.18*1e-2;
delta_Ns = 0.8*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns) .* sqrt(s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf);
G_Nv = double(G_Sv);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
%the integral for H_Sv is going to make it into the final expression. However, matlabFunction cannot work
%with vpaintegral() without falling over, and int() takes a long time to figure out that it is not going to
%be able to find an answer. In the interests of performance, code vpaintegral() here, and we will use
%a hack to swap it out to int() later.
H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf);   %this takes a long time to figure out it cannot do anything.
%%lambda_Ss= c_Ns./fr_Nv;
p_Nv = c_Ns ./ fr_Nv;
%N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, lambda_Ss);    %Cannot substitute and replace the lambda by o and I don't know why
%%Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, lambda_Ss);
%{
G_Sv / G_Nv are independent of any variables so no point doing the subs()
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Nv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
%}
N_Nv = G_Nv;
%{
H_Sv is independent of lambda, so there is no point doing the subs()
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %Same substitution problem
%}
Q_Sv = H_Sv;
temp_Sv = (1/(4 .* mur_Ns)) .* (N_Nv ./ Q_Sv);
SEHth_Sv = sqrt(real(temp_Sv).^2 + imag(temp_Sv).^2);   %avoid using abs, it is not continuously differentiable
% %% Se
Se_Nv = [12, 19, 21, 26, 29 ,34 , 38 ,39 ,41, 42 ,41, 42];
%{
output_Sv is not used so no point computing it
output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);
%}
%X= abs((Se_Nv - SEHth_Sv).^2);     %abs is not differentiable, use the smooth function instead
%{
SEHth_Sv is already abs() so we do not have to worry about imaginary
%X = (real(Se_Nv - SEHth_Sv)).^2+(imag(Se_Nv - SEHth_Sv)).^2;
%}
X = (Se_Nv - SEHth_Sv).^2;
%{
X is independent of lambda, so no point doing the subs
F=sum(subs(X,lambda_Ss,p_Nv));
%}
F = sum(X);
Fint = mapSymType(F, 'vpaintegral', @(V) feval(symengine, 'hold(int)', feval(symengine, 'op', V, 1), 'lambda = 0..infinity'));
%% objective harmonization
Fc = matlabFunction(Fint,'Vars',{s_Ss});   %if you send F to file then optimize must be off
gama_fun = matlabFunction(gama_Sv, 'Vars', {s_Ss}, 'File', 'gama_fun.m', 'optimize', true);
obj = Fc;
A=[]; 
b=[];
s_Ss0=1e-22;
lb=1e-22;
ub=1e8;
Aeq = [];
beq = [];
% %nonlcon = @(s_Ss) constrain(s_Ss, fr_Nv, muo_Ns, mur_Ns, gama_fun);
nonlcon= @(s_Ss) constrain(s_Ss);
% optimize with fmincon
%[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN]
%%fmin= fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS);
%%fmin = fmincon(obj, s_Ss0, A, b, Aeq, beq, lb, ub, nonlcon);
s_Ss = fmincon(obj,s_Ss0,A,b,Aeq,beq,lb,ub,nonlcon);

% show final objective
disp(['Final Objective: ' num2str(obj(s_Ss))])

%Constrain function
function [c,ceq] = constrain(s_Ss)
c = [];
ceq = [];
end

[Walter Roberson]

It looks to me as if the minimum is always at your lower bound.

[Salwa Ben Mbarek]

Thank you sir, I see.

I've changed the lower and upper bound to -inf, +inf and the value of the starting point to see if this has impact on the minimum. I have the same result : 17427189.0848 in the program below . Is it mathematically wrong or it can happen that the minimum is at the lower bound for a program like this? I mean is it problematic if the minimum is always at the lower bound?

Thank you very much !

The program is :

% Sth
z_Ns = 92.18*1e-2;
delta_Ns = 0.8*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns) .* sqrt(s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf);
G_Nv = double(G_Sv);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
%the integral for H_Sv is going to make it into the final expression. However, matlabFunction cannot work
%with vpaintegral() without falling over, and int() takes a long time to figure out that it is not going to
%be able to find an answer. In the interests of performance, code vpaintegral() here, and we will use
%a hack to swap it out to int() later.
H_Sv = vpaintegral(M_Sv, lambda_Ss , 0, inf);   %this takes a long time to figure out it cannot do anything.
%%lambda_Ss= c_Ns./fr_Nv;
p_Nv = c_Ns ./ fr_Nv;
%N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, lambda_Ss);    %Cannot substitute and replace the lambda by o and I don't know why
%%Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, lambda_Ss);
%{
G_Sv / G_Nv are independent of any variables so no point doing the subs()
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Nv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
%}
N_Nv = G_Nv;
%{
H_Sv is independent of lambda, so there is no point doing the subs()
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %Same substitution problem
%}
Q_Sv = H_Sv;
temp_Sv = (1/(4 .* mur_Ns)) .* (N_Nv ./ Q_Sv);
SEHth_Sv = sqrt(real(temp_Sv).^2 + imag(temp_Sv).^2);   %avoid using abs, it is not continuously differentiable
% %% Se
Se_Nv = [12, 19, 21, 26, 29 ,34 , 38 ,39 ,41, 42 ,41, 42];
%{
output_Sv is not used so no point computing it
output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);
%}
%X= abs((Se_Nv - SEHth_Sv).^2);     %abs is not differentiable, use the smooth function instead
%{
SEHth_Sv is already abs() so we do not have to worry about imaginary
%X = (real(Se_Nv - SEHth_Sv)).^2+(imag(Se_Nv - SEHth_Sv)).^2;
%}
X = (Se_Nv - SEHth_Sv).^2;
%{
X is independent of lambda, so no point doing the subs
F=sum(subs(X,lambda_Ss,p_Nv));
%}
F = sum(X);
%%Fint = mapSymType(F, 'vpaintegral', @(V) feval(symengine, 'hold(int)', feval(symengine, 'op', V, 1), 'lambda= 3.0000..0.0010'));
Fint = mapSymType(F, 'vpaintegral', @(V) feval(symengine, 'hold(int)', feval(symengine, 'op', V, 1), 'lambda = 0..infinity'));
%% objective harmonization
Fc = matlabFunction(Fint,'Vars',{s_Ss});   %if you send F to file then optimize must be off
gama_fun = matlabFunction(gama_Sv, 'Vars', {s_Ss}, 'File', 'gama_fun.m', 'optimize', true);
obj = Fc;
A=[];
b=[];
% s_Ss0=1e-22;
% lb=1e8;
% ub=1e-22;
s_Ss0=-2;
lb=-inf;
ub=inf;
%ub=1e8;
Aeq = [];
beq = [];
% %nonlcon = @(s_Ss) constrain(s_Ss, fr_Nv, muo_Ns, mur_Ns, gama_fun);
nonlcon= @(s_Ss) constrain(s_Ss);
% optimize with fmincon
%[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN]
%%fmin= fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS);
%%fmin = fmincon(obj, s_Ss0, A, b, Aeq, beq, lb, ub, nonlcon);
s_Ss = fmincon(obj,s_Ss0,A,b,Aeq,beq,lb,ub,nonlcon);

% show final objective
disp(['Final Objective: ' num2str(obj(s_Ss))])

Constrain function :

function [c,ceq] = constrain(s_Ss)
c = [];
ceq = [];
end

[Walter Roberson]

I switched over to symbolic and have been trying to find a more precise value for 1e-20, in order to get an idea of how much double precision is limiting the solution. The calculation has been going for quite a number of hours so far...

[Salwa Ben Mbarek]

Thank you so much sir for your precious help. What variable did you put on symbolic ? Do not bother yourself if the calculation takes all of this Time .. it's already very kind of you to help me ..I Can try to run it on my laptop don't worry . Sir ,did you use the same program to find a more precise value for 1e-20 ? Is the double precision generated automatically or we should introduce the syntax "double" before the calculation of the minimum ? Thank you very much !

[Walter Roberson]

If you take Fint and change the vpaintegral() to int(), and then try to evaluate the resulting Fint at sym(10)^(-20) to get a symbolic solution, then it takes a long long time.

It might be worth testing with the Fint that uses vpaintegral() like we already have, but tell vpaintegral() to use RelTol of about 1e-30 or so.

The hypothesis behind these tests would be that perhaps there are substantial differences in output values for inputs less than 1e-9, but that the conversion to double precision done by matlabFunction, and the errors introduced by order of operations, are masking the true responses.

However, in all of the tests I have done so far, it has seemed pretty consistent that once you get beyond a threshold of about 1e-9 where all of the answers come out the same to the limits of double precision, that the output rises fairly rapidly.

I have not tried evaluating over larger ranges (say up to 100) to see whether perhaps there is an oscillation .

[Salwa Ben Mbarek]

Hello Sir,

Thank you so much for your help . I understand.

I told vpainintegral to use Reltol , this is the program :

% Sth
z_Ns = 92.18*1e-2;
delta_Ns = 0.8*1e-3;
c_Ns = 3e8;
muo_Ns = 4 .* pi .* 10^-7;
fr_Nv = [10e3, 20e3, 50e3, 100e3, 200e3, 500e3, 1e6, 2e6, 5.1e6, 10.1e6, 20.1e6, 30.1e6];
mur_Ns = 1;
w_Nv = 2 .* pi .* fr_Nv;
%in this section, I chose to put lambda as a sym in calculations, and to
%substitute its value at the end
lambda_Ss = sym('lambda');
gamao_Nv = 1j.*w_Nv./c_Ns;
%s is symbolic
s_Ss = sym('s');
gama_Sv = sqrt(1j .* 2 .* pi .* fr_Nv .* muo_Ns .* mur_Ns) .* sqrt(s_Ss);  %uses s
too_Sv = sqrt(lambda_Ss.^2 - gamao_Nv.^2);                           %does not use s
to_Sv = sqrt(lambda_Ss.^2 - gama_Sv.^2);                             %uses s
K_Sv = ((to_Sv./too_Sv + mur_Ns).^2 - (to_Sv./too_Sv - mur_Ns).^2 .* exp(-2 .* to_Sv .* delta_Ns)).^-1; %uses s
R_Ns = 16.05*1e-2;                                                   %does not use s
J_Sv = besselj(1, lambda_Ss .* R_Ns);                                %does not use s
L_Sv = lambda_Ss .* lambda_Ss .* ((too_Sv).^-1) .* J_Sv .* exp(-too_Sv .* z_Ns);   %does not use s
% L is now symbolic
%now notice that here we are integrating the purely numeric L with respect to symbolic lambda
G_Sv = vpaintegral(L_Sv, lambda_Ss, 0, inf,'RelTol',1e-30,'AbsTol', 0);
G_Nv = double(G_Sv);
M_Sv = K_Sv .* (lambda_Ss.^2) .* (to_Sv) .* (too_Sv.^-2) .* J_Sv .* exp(-too_Sv .* z_Ns - (to_Sv - too_Sv)) .* delta_Ns;
%the integral for H_Sv is going to make it into the final expression. However, matlabFunction cannot work
%with vpaintegral() without falling over, and int() takes a long time to figure out that it is not going to
%be able to find an answer. In the interests of performance, code vpaintegral() here, and we will use
%a hack to swap it out to int() later.
H_Sv = vpaintegral(M_Sv, lambda_Ss, 0, inf, 'RelTol',1e-30,'AbsTol', 0);   %this takes a long time to figure out it cannot do anything.
%%lambda_Ss= c_Ns./fr_Nv;
p_Nv = c_Ns ./ fr_Nv;
%N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Sv, lambda_Ss);    %Cannot substitute and replace the lambda by o and I don't know why
%%Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, lambda_Ss);
%{
G_Sv / G_Nv are independent of any variables so no point doing the subs()
N_Sv = arrayfun(@(G,P) subs(G, lambda_Ss, P), G_Nv, p_Nv);    %Cannot substitute and replace the lambda by o and I don't know why
%}
N_Nv = G_Nv;
%{
H_Sv is independent of lambda, so there is no point doing the subs()
Q_Sv = arrayfun(@(H,P) subs(H, lambda_Ss, P), H_Sv, p_Nv);    %Same substitution problem
%}
Q_Sv = H_Sv;
temp_Sv = (1/(4 .* mur_Ns)) .* (N_Nv ./ Q_Sv);
SEHth_Sv = sqrt(real(temp_Sv).^2 + imag(temp_Sv).^2);   %avoid using abs, it is not continuously differentiable
% %% Se
Se_Nv = [12, 19, 21, 26, 29 ,34 , 38 ,39 ,41, 42 ,41, 42];
%{
output_Sv is not used so no point computing it
output_Sv = arrayfun(@(S, SE) subs(S, s_Ss, SE), SEHth_Sv, Se_Nv);
%}
%X= abs((Se_Nv - SEHth_Sv).^2);     %abs is not differentiable, use the smooth function instead
%{
SEHth_Sv is already abs() so we do not have to worry about imaginary
%X = (real(Se_Nv - SEHth_Sv)).^2+(imag(Se_Nv - SEHth_Sv)).^2;
%}
X = (Se_Nv - SEHth_Sv).^2;
%{
X is independent of lambda, so no point doing the subs
F=sum(subs(X,lambda_Ss,p_Nv));
%}
F = sum(X);
%%Fint = mapSymType(F, 'vpaintegral', @(V) feval(symengine, 'hold(int)', feval(symengine, 'op', V, 1), 'lambda= 3.0000..0.0010'));
Fint = mapSymType(F, 'vpaintegral', @(V) feval(symengine, 'hold(int)', feval(symengine, 'op', V, 1), 'lambda = 0..infinity'));

%%RelTol',1e-30;
%%% objective harmonization
Fc = matlabFunction(Fint,'Vars',{s_Ss});   %if you send F to file then optimize must be off
gama_fun = matlabFunction(gama_Sv, 'Vars', {s_Ss}, 'File', 'gama_fun.m', 'optimize', true);
obj = Fc;
A=[];
b=[];
% s_Ss0=1e-22;
% lb=1e8;
% ub=1e-22;
s_Ss0=-2;
lb=-inf;
ub=inf;
%ub=1e8;
Aeq = [];
beq = [];
% %nonlcon = @(s_Ss) constrain(s_Ss, fr_Nv, muo_Ns, mur_Ns, gama_fun);
nonlcon= @(s_Ss) constrain(s_Ss);
% optimize with fmincon
%[X,FVAL,EXITFLAG,OUTPUT,LAMBDA,GRAD,HESSIAN]
%%fmin= fmincon(FUN,X0,A,B,Aeq,Beq,LB,UB,NONLCON,OPTIONS);
%%fmin = fmincon(obj, s_Ss0, A, b, Aeq, beq, lb, ub, nonlcon);
s_Ss = fmincon(obj,s_Ss0,A,b,Aeq,beq,lb,ub,nonlcon);
% show final objective
disp(['Final Objective: ' num2str(obj(s_Ss))])

Constrain function :

function [c,ceq] = constrain(s_Ss)
c = [];
ceq = [];
end

I still get the value of 1.7427e+07. You're right, the value increases rapidly

obj(100) =3.6222e+94
obj(1000) = 6.9310e+297

Do you think is it logical to have this minimum even when we use Reltol and have -inf and + inf as lower and upper bounds?

Thank you so much !

[Walter Roberson]

Whether it makes sense or not, with the equations you give, the derivative of the function is nearly linear up to about 2E-8, falling a little less than linear after that. The derivative is approximately 4E10*s for values between 0 and 2E-8. In the tracking I did, the derivative for positive s is always positive.

... And that means that the minimum value of the function is going to be at the lowest point that you evaluate the function at.