The conflict between 'ArrayValued' option in integration and specia... (2024)

Table of Contents
Answers (0) See Also Categories Tags Products Release Community Treasure Hunt

2 views (last 30 days)

Show older comments

tsukatsuki on 10 Apr 2024

  • Link

    Direct link to this question

    https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x

⋮

  • Link

    Direct link to this question

    https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x

Edited: Torsten on 11 Apr 2024

Open in MATLAB Online

I wrote a function HCDd to batch process some integrals and put them in a matrix, but encountered problems when it came to the special function psi (X) (digamma function).

When the 'ArrayValue' option is set to True, the 'integral' function seems to only accept input of single and double floating-point types. But the integration process involves negative input to the function psi (X) (see my custom function Rnd), and in the case of non-symbolic input, psi (X) will call the built-in function of Matlab instead of the function of the same name in the Symbolic Math Toolbox. But the built-in function psi (X) can only accept non-negative real input, which breaks my calculations.

I want to optimize computational efficiency as much as possible, so I am not inclined to use Arrayfun or for loops. Is there a better solution?

digits(10);

HCDd(4, 0, 1) %for test

function omega_t = omega(t)

omega_t = vpa(3 + ((t - 2)^4 - 8 * (t - 2)^2) / 8);

end

function eg = Eg(n, a0)

efun = @(egtmp) sqrt(2) * gamma(-egtmp/2 + 3/4) / gamma(-egtmp/2 + 1/4) - 1/a0;

if n~= -1

if a0 ~= 0

smin = 2 * n + 0.5;

smax = 2 * n + 2.5;

else

eg = vpa(1.5 + 2 * n);

return;

end

else

if a0 ~= 0

syms egtmp1;

efun1 = sqrt(2) * gamma(-egtmp1/2 + 3/4) / gamma(-egtmp1/2 + 1/4) - 1/a0 == 0;

eg = vpasolve(efun1,-1);

return;

else

eg = -inf;

return;

end

end

eg = vpa(GlobalSolve(efun,1,smin,smax)); %GlobalSolve is a custom function unrelated to this problem, used for finding roots

end

function nu0 = nu(n, a0)

nu0 = Eg(n, a0) / 2 - 3/4;

end

function rnd = Rnd(nu_n_a0, r, t) %here use the fuction psi(X)

rnd = omega(t)^(3/4) .* 2 .^ (nu_n_a0 + 3/2) .* pi .^ (-1/4) ...

.* sqrt(gamma(-2 .* nu_n_a0 - 1) ./ (psi(-nu_n_a0) - psi(-nu_n_a0 - 1/2))) ...

.* exp(-omega(t) / 2 .* r .^ 2) ...

.* kummerU(-nu_n_a0, 3/2, omega(t) .* r .^ 2);

end

function nuarray = nuArray(M, a0)

nuarray = vpa(zeros(1, M));

for m = 1:M

nuarray(1, m) = nu(m-2, a0);

end

end

function hcd = HCDd(M, t, a0) %problem in here

syms nutmp rtmp r nu_mAtmp nu_nAtmp

nu_A = nuArray(M, a0);

nu_mA = double(repmat(nu_A, M, 1));

nu_nA = nu_mA.';

dRnd = matlabFunction(diff(Rnd(nutmp, rtmp, t), rtmp));

func = Rnd(nu_mAtmp, r, t) .* 2 .* r .^3 .* dRnd(nu_nAtmp, r) ...

+ Rnd(nu_mAtmp, r, t) .* 3 .* r .^2 .* Rnd(nu_nAtmp, r, t);

fun = matlabFunction(func, 'vars', {r nu_mAtmp nu_nAtmp} );

hcd = integral(@(r) fun(r, nu_mA, nu_nA), 0, inf, 'ArrayValued', true);

end

6 Comments

Show 4 older commentsHide 4 older comments

Torsten on 10 Apr 2024

Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127256

See Also
Is there a way to vpasolve elementwise?Using fzero with syms or finding an alternative to syms

⋮

  • Link

    Direct link to this comment

    https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127256

If you set "ArrayValued",true in the call to integral, only one numerical value for r is passed to "fun" once at a time. So this has nothing to do with symbolic or numerical evaluation of the "psi" function (as far as I can see).

But in order to help, it would be easiest to include executable code, i.e. supplying "GlobalSolve" and other functions that might be necessary.

tsukatsuki on 10 Apr 2024

Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127411

⋮

  • Link

    Direct link to this comment

    https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127411

Edited: tsukatsuki on 10 Apr 2024

Open in MATLAB Online

@Torsten Please note that the parameters of the integration process come from the two matrices nu_mA and nu_nA, and r is the integration variable here. I have sorted out my two different attempts on the HCDd function and the corresponding error messages, and also attached the content of the GlobalSolve function you mentioned.

First I run the original program above, Matlab says

Error using psi

X must be nonnegative.

Error in

symengine>@(r,nu_mAtmp,nu_nAtmp)2.0.^(nu_mAtmp+3.0./2.0).*r.^3.*exp(r.^2.*-5.0e-1).*(2.0.^(nu_nAtmp+3.0./2.0).*r.*exp(r.^2.*(-1.0./2.0)).*((1.0./r.^2.*kummerU(-nu_nAtmp,3.0./2.0,r.^2))./2.0-1.0./r.^2.*(nu_nAtmp+1.0./2.0).*kummerU(-nu_nAtmp,1.0./2.0,r.^2)).*sqrt(gamma(nu_nAtmp.*-2.0-1.0)./(psi(-nu_nAtmp)-psi(-nu_nAtmp-1.0./2.0))).*1.502251088929885+2.0.^(nu_nAtmp+3.0./2.0).*r.*exp(r.^2.*(-1.0./2.0)).*sqrt(gamma(nu_nAtmp.*-2.0-1.0)./(psi(-nu_nAtmp)-psi(-nu_nAtmp-1.0./2.0))).*kummerU(-nu_nAtmp,3.0./2.0,r.^2).*7.511255444649425e-1).*sqrt(gamma(nu_mAtmp.*-2.0-1.0)./(psi(-nu_mAtmp)-psi(-nu_mAtmp-1.0./2.0))).*kummerU(-nu_mAtmp,3.0./2.0,r.^2.*1.0).*-1.502251088929885+2.0.^(nu_mAtmp+3.0./2.0).*2.0.^(nu_nAtmp+3.0./2.0).*r.^2.*exp(r.^2.*-1.0).*sqrt(gamma(nu_mAtmp.*-2.0-1.0)./(psi(-nu_mAtmp)-psi(-nu_mAtmp-1.0./2.0))).*sqrt(gamma(nu_nAtmp.*-2.0-1.0)./(psi(-nu_nAtmp)-psi(-nu_nAtmp-1.0./2.0))).*kummerU(-nu_mAtmp,3.0./2.0,r.^2.*1.0).*kummerU(-nu_nAtmp,3.0./2.0,r.^2.*1.0).*1.692568750643269

Error in Copy_of_zero3d>@(r)fun(r,nu_mA,nu_nA) (line 61)

hcd = vpa(integral(@(r) fun(r, nu_mA, nu_nA), 0, inf, 'ArrayValued', true));

Error in integralCalc/iterateArrayValued (line 156)

fxj = FUN(t(1)).*w(1);

Error in integralCalc/vadapt (line 130)

[q,errbnd] = iterateArrayValued(u,tinterval,pathlen);

Error in integralCalc (line 83)

[q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral (line 87)

Q = integralCalc(fun,a,b,opstruct);

Error in Copy_of_zero3d>HCDd (line 61)

hcd = vpa(integral(@(r) fun(r, nu_mA, nu_nA), 0, inf, 'ArrayValued', true));

Error in Copy_of_zero3d (line 4)

HCDd(4, 0, 1)

Then I remove the double function on nu_mA and nu_nA to restore them to sym type, like this

function hcd = HCDd(M, t, a0)

syms nutmp rtmp r nu_mAtmp nu_nAtmp

nu_A = nuArray(M, a0);

nu_mA = repmat(nu_A, M, 1);

nu_nA = nu_mA.';

dRnd = matlabFunction(diff(Rnd(nutmp, rtmp, t), rtmp))

func = Rnd(nu_mAtmp, r, t) .* 2 .* r .^3 .* dRnd(nu_nAtmp, r) ...

+ Rnd(nu_mAtmp, r, t) .* 3 .* r .^2 .* Rnd(nu_nAtmp, r, t)

fun = matlabFunction(func, 'vars', {r nu_mAtmp nu_nAtmp} )

hcd = vpa(integral(@(r) fun(r, nu_mA, nu_nA), 0, inf, 'ArrayValued', true));

end

then Matlab says

Error using superiorfloat

Inputs must be floats, namely single or double.

Error in integralCalc/iterateArrayValued (line 159)

outcls = superiorfloat(x,fxj);

Error in integralCalc/vadapt (line 130)

[q,errbnd] = iterateArrayValued(u,tinterval,pathlen);

Error in integralCalc (line 83)

[q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral (line 87)

Q = integralCalc(fun,a,b,opstruct);

Error in Copy_of_zero3d>HCDd (line 61)

hcd = vpa(integral(@(r) fun(r, nu_mA, nu_nA), 0, inf, 'ArrayValued', true));

Error in Copy_of_zero3d (line 4)

HCDd(4, 0, 1)

So I think the problem occurs as I mentioned at the beginning.

Alternatively, here's what GlobalSolve.m is all about, and you can now try running the program in the question on your own Matlab.

function varargout = GlobalSolve(varargin)

n = varargin{2};

obj_fun = varargin{1};

start_num = 100;

if(nargin==2)

lb = -inf(1,n);

ub = inf(1,n);

end

if(nargin==4)

lb = varargin{3};

ub = varargin{4};

end

if(nargin==5)

if isempty(varargin{3})

lb = -inf(1,n);

else

lb = varargin{3};

end

if isempty(varargin{4})

ub = inf(1,n);

else

ub = varargin{4};

end

start_num = varargin{5};

end

options = optimoptions('particleswarm','display','off');

x0 = particleswarm(@(x)sum(obj_fun(x).^2),n,lb,ub,options);

opts = optimoptions(@fmincon,'Algorithm','interior-point','Display','off');

problem = createOptimProblem('fmincon','x0',x0,'objective',@(x)0,'lb',lb,'ub',ub,'nonlcon',@(x)constrain(x,obj_fun));

gs = MultiStart('FunctionTolerance',1e-2,'XTolerance',1e-2);

[x,~,exitflag,~,allmins] = gs.run(problem,start_num);

fval = obj_fun(x);

if n==1

for ii = 1:numel(allmins)

multiSol(ii) = allmins(ii).X;

end

multiSol = sort(multiSol);

for ii = 1:numel(allmins)

multiFval{ii} = obj_fun(multiSol(ii));

end

else

for ii = 1:numel(allmins)

multiSol{ii} = allmins(ii).X;

multiFval{ii} = obj_fun(multiSol{ii});

end

end

switch nargout

case 1

varargout{1} = x;

case 2

varargout{1} = x;

varargout{2} = fval;

case 3

varargout{1} = x;

varargout{2} = fval;

varargout{3} = exitflag;

case 4

varargout{1} = x;

varargout{2} = fval;

varargout{3} = exitflag;

varargout{4} = multiSol;

case 5

varargout{1} = x;

varargout{2} = fval;

varargout{3} = exitflag;

varargout{4} = multiSol;

varargout{5} = multiFval;

end

function [c ceq] = constrain(x,f)

ceq = f(x);

c =[];

tsukatsuki on 10 Apr 2024

Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127421

⋮

  • Link

    Direct link to this comment

    https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127421

Open in MATLAB Online

Also, I tried a less efficient one using a for loop. In fact, the following version of the HCDd function can run well, but this does not solve the optimization problem.

function hcd = HCDd(M, t, a0)

syms nutmp rtmp r

nu_A = nuArray(M, a0);

hcd = vpa(zeros(M));

dRnd = matlabFunction(diff(Rnd(nutmp, rtmp, t), rtmp))

for m = 1:M

for n = 1:M

func = Rnd(nu_A(1, m), r, t) .* 2 .* r .^3 .* dRnd(nu_A(1, n), r) ...

+ Rnd(nu_A(1, m), r, t) .* 3 .* r .^2 .* Rnd(nu_A(1, n), r, t)

fun = matlabFunction(func, 'vars', {r} )

hcd(m, n) = vpa(integral(fun, 0, inf));

end

end

end

Torsten on 10 Apr 2024

Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127741

⋮

  • Link

    Direct link to this comment

    https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127741

Edited: Torsten on 10 Apr 2024

Also, I tried a less efficient one using a for loop. In fact, the following version of the HCDd function can run well, but this does not solve the optimization problem.

What do you mean by "does not solve the optimization problem" ?

Do you mean that the code is not optimized or that the solution you get is wrong in your opinion ?

tsukatsuki on 11 Apr 2024

Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127981

⋮

  • Link

    Direct link to this comment

    https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3127981

@Torsten Sorry I used the wrong words, I'm not a native English speaker. What I'm trying to say is that using a for loop slows down the running of the function, which in the final program may be run thousands of times with a larger matrix order M. So I want to use a vectorization way.

Torsten on 11 Apr 2024

Direct link to this comment

https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3128296

⋮

  • Link

    Direct link to this comment

    https://matlabcentral.mathworks.com/matlabcentral/answers/2105101-the-conflict-between-arrayvalued-option-in-integration-and-special-function-psi-x#comment_3128296

Edited: Torsten on 11 Apr 2024

Open in MATLAB Online

I'd generate the function "fun" once before your "big simulation", write it to file and replace calls to psi(x) by double(psi(vpa(x))) manually.

x = -0.625;

double(psi(vpa(x)))

ans = -1.1540

psi(x)

Error using psi
X must be nonnegative.

I don't know if the symbolic psi-function will be drastically slower than the numerical one.

Sign in to comment.

Sign in to answer this question.

Answers (0)

Sign in to answer this question.

See Also

Categories

Mathematics and OptimizationSymbolic Math ToolboxMathematicsCalculus

Find more on Calculus in Help Center and File Exchange

Tags

  • special function
  • psi
  • integration
  • numerical integration
  • symbolic

Products

  • MATLAB
  • Symbolic Math Toolbox

Release

R2023b

Community Treasure Hunt

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

Start Hunting!

An Error Occurred

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.


The conflict between 'ArrayValued' option in integration and specia... (8)

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

Americas

  • América Latina (Español)
  • Canada (English)
  • United States (English)

Europe

  • Belgium (English)
  • Denmark (English)
  • Deutschland (Deutsch)
  • España (Español)
  • Finland (English)
  • France (Français)
  • Ireland (English)
  • Italia (Italiano)
  • Luxembourg (English)
  • Netherlands (English)
  • Norway (English)
  • Österreich (Deutsch)
  • Portugal (English)
  • Sweden (English)
  • Switzerland
    • Deutsch
    • English
    • Français
  • United Kingdom(English)

Asia Pacific

Contact your local office

The conflict between 'ArrayValued' option in integration and specia... (2024)
Top Articles
The Post-Crescent from Appleton, Wisconsin
naturalisme, Algemeen letterkundig lexicon - DBNL
2020 - UPDATED: May 2024
Southern Caribbean from Baltimore, Royal Caribbean, 29th March 2025 – Planet Cruise
Connect America Fund, Connect America Fund Phase II Auction, The Uniendo a Puerto Rico Fund and the Connect USVI Fund, Rural Digital Opportunity Fund, Rural Digital Opportunity Fund Auction, Establishing a 5G Fund for Rural America, Letters of Credit for
Getting Started - Emergency Connectivity Fund - USAC.ORG
Milli Vanilli: Where Are The Band Members Now?
The Tragic Death Of Rob Pilatus From Milli Vanilli - Nicki Swift
NYT Mini Crossword: 8 Strategies To Solve It Fast
Renovo Record Live Stream
Pacemaker | MedUni Vienna
Frontiers | Effects of Rhythm and Rate-Controlling Drugs in Patients With Permanent His-Bundle Pacing
Latest Posts
ABC News – Breaking News, Latest News and Videos
Best CBD Gummies for March 2024
Article information

Author: Manual Maggio

Last Updated:

Views: 5691

Rating: 4.9 / 5 (49 voted)

Reviews: 88% of readers found this page helpful

Author information

Name: Manual Maggio

Birthday: 1998-01-20

Address: 359 Kelvin Stream, Lake Eldonview, MT 33517-1242

Phone: +577037762465

Job: Product Hospitality Supervisor

Hobby: Gardening, Web surfing, Video gaming, Amateur radio, Flag Football, Reading, Table tennis

Introduction: My name is Manual Maggio, I am a thankful, tender, adventurous, delightful, fantastic, proud, graceful person who loves writing and wants to share my knowledge and understanding with you.