construct a local function using a mixture of real and symbolic math
1 view (last 30 days)
Show older comments
Would love to have help with two things. This local function runs, but …
- I’m getting partially symbolic output and I need real-valued output. (See output sample at bottom.)
- How to make this more efficient. This version runs in about 20 seconds, but I really want to set finer grained values in the loops; i.e., generating 100 times the current number of function calls and that will burn some time in the present form. Is it possible to turn symsum into a @ type function. I tried and failed.
Thanks so much if you can help on this. This is my first venture with the symbolic toolbox, so my guess is that I am really screwing this up and missing important pieces.
Dave C.
Here’s the function “symbolically”; G = gamma function
Sum(kk=0 to Inf) { (rho*sqrt(ab)/(1-rho^2) * [G((kk+1)/2]/[kk! G((kk+m)/2)] }
% Contours_ChiSq.m
% In progress
% ISum is the local function that needs work.
clear all; close all; clc;
rho = 0.5;
k = 1;
for a = 0.1:0.05:1
for b = (1-a):0.05:1
Xab(k,:) = [a b ISum(rho,a,b)];
k = k+1;
end
end
Xab;
%--------------------------------------------------------------------------
% Call the LOCAL infinite sum function below
%--------------------------------------------------------------------------
sympref('FloatingPointOutput',true);
function [s] = ISum(rho,a,b)
m = 4;
kk =sym('kk');
S = (rho*sqrt(a*b)/(1-rho))^kk;
Ng = gamma((kk+1)/2);
Dg = factorial(kk)*gamma((kk+m)/2);
ss = vpa(symsum(S,kk,0,inf),5);
s = ss*vpa(Ng/Dg,5);
end
Sample output:
Xab =
[0.1000, 0.9000, (1.4286*gamma(0.5000*kk + 0.5000))/(gamma(0.5000*kk + 2)*factorial(kk))]
[0.1000, 0.9500, (1.4455*gamma(0.5000*kk + 0.5000))/(gamma(0.5000*kk + 2)*factorial(kk))]
[0.1000, 1, (1.4625*gamma(0.5000*kk + 0.5000))/(gamma(0.5000*kk + 2)*factorial(kk))]
[0.1500, 0.8500, (1.5554*gamma(0.5000*kk + 0.5000))/(gamma(0.5000*kk + 2)*factorial(kk))]
5 Comments
Torsten
on 26 Dec 2023
Edited: Torsten
on 26 Dec 2023
S = (rho*sqrt(a*b)/(1-rho^2))^kk;
instead of
S = (rho*sqrt(a*b)/(1-rho))^kk;
and
ss = vpa(symsum(S*Ng/Dg,kk,0,inf),5);
instead of
ss = vpa(symsum(S,kk,0,inf),5);
(Ng and Dg depend on kk - they cannot be taken out of the sum).
Why do you work with symbolic variables at all ? Sum with doubles until the size of the summands get smaller than a specified tolerance.
Answers (1)
Walter Roberson
on 26 Dec 2023
Edited: Walter Roberson
on 26 Dec 2023
You can simplify a lot.
rho = sym(0.5);
syms a b
assume(a > 0 & a <= 1);
assumeAlso(b > 0 & b <= 1);
m = 4;
kk =sym('kk');
S = (rho*sqrt(a*b)/(1-rho))^kk;
Ng = gamma((kk+1)/2);
Dg = factorial(kk)*gamma((kk+m)/2);
SS = symsum(S,kk,0,inf)
The second condition will not occur.
So your ss will be infinity if a = 1 and b = 1, and will be 1/(1 - sqrt(a)*sqrt(b)) otherwise. For practical purposes, you can reduce that to just 1/(1 - sqrt(a)*sqrt(b)) as that will come out as infinity when a = 1 and b = 1.
See Also
Categories
Find more on Assumptions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!