Symsum function with odd number and infinity upper limit

Hi,

I'm trying to do a sum series of the following equation with odd numbers and infinity upper limit.

When running the following code, i have a littteral value not a unique value which display in the command window. How can i do to have a unique value in a variable? How can i define the calculation with infinity upper limit and odd number? Please

The aim is to be able to plot the distribution of p over x and y.

%% Analytical calculation

clear all
close all
clc


%% Input data

f = 15; % frequency [Hz]

omega = 2*pi.*f; % angular frequency

ho = 0.3*1e-3/2; % air gap height [m]

mu = 18.1e-6; % dynamic viscosity [Ns/m2]

rhoo = 1.206; % mean density [Kg/m3]

Co = 343; % undisturbed speed of sound [m/s]

Cp = 1004; % specific heat at constant pressure [J/KgK]

Cv = 716; % specific heat at constant volume [J/KgK]

lambda = 25.6e-3; % thermal conductivity [W/mK]

lx = 250*1e-3/2; % half length of plate [m]

ly = 180*1e-3/2; % half width of plate [m]

h = 0.05.*1e-3; % displacement amplitude

po = 1.01e5; % ambient pressure [N/m2]

%% Dimensionless Parameters

s = ho.*sqrt(rhoo.*omega./mu); % shear wave number

k = omega.*ho./Co; % reduced frequency

sigma = sqrt(mu.*Cp./lambda); % square root of Prandtl number

g = ho./lx; % nawroness of the gap

gamma = Cp/Cv; % ratio of specific heats

a = ly/lx; % aspect ratio of the panel

kx = omega*lx/Co; % normalised wave number along x-axis

ky = omega*ly/Co; % normalised wave number along y-axis


Bs =(tanh(s.*sqrt(i))./(s.*sqrt(i)))-1;

Bssigma =(tanh(s*sigma.*sqrt(i))./(s*sigma.*sqrt(i)))-1;

nssigma = (1+((gamma-1)/gamma).*Bssigma).^(-1);

Gamma = sqrt(gamma/(nssigma*Bs));


Ceffs = Co.*sqrt(nssigma.*Bs/gamma); % effective speed sound [m/s]


syms x y q 

x = 1

x = 1

y = 1

y = 1

D = sqrt(((q*pi/(2*a))^2) - (kx*Gamma)^2);

SUMp=0;
a=1:2:9;
 for jj=1:length(a)
     
     p = (2*nssigma*Gamma^2*kx*ky*h/pi).*symsum((((-1)^((q-1)/2)/(q*(D)^2)))*((cosh(D*x/kx)/cosh(D))-1)*cos(q*pi*y/(2*ky)),q,a(jj),a(jj));     

    SUMp=SUMp+p;
end
pretty(SUMp);

   / 1715 \ /                / 1372 #5 \ \ /      114581679025877            380147757317807      \
cos| ---- | | cosh(#5) - cosh| ------- | | | ------------------------ - ------------------------i |
   \   3  / \                \  15 pi  / / \ 885443715538058477568000   166020696663385964544000  /
---------------------------------------------------------------------------------------------------
                             /   2      9500033100991       253138413248751   \
                    cosh(#5) | pi  + ------------------ - ------------------i |
                             \       225179981368524800   281474976710656000  /

        / 1715 \ /                / 1372 #4 \ \ /       1031235111232893           1140443271953421     \
     cos| ---- | | cosh(#4) - cosh| ------- | | | - ------------------------ + -----------------------i |
        \   9  / \                \  15 pi  / / \   295147905179352825856000   18446744073709551616000  /
   + ----------------------------------------------------------------------------------------------------
                                   /   2     85500297908919      2278245719238759   \
                          cosh(#4) | pi  + ------------------ - ------------------i |
                                   \       225179981368524800   281474976710656000  /

        / 1715 \ /                / 1372 #3 \ \ /     27843348003288111         30791968342742367     \
     cos| ---- | | cosh(#3) - cosh| ------- | | | ------------------------ - -----------------------i |
        \  27  / \                \  15 pi  / / \ 295147905179352825856000   18446744073709551616000  /
   + --------------------------------------------------------------------------------------------------
                                  /   2     769502681180271     20504211473148831  \
                         cosh(#3) | pi  + ------------------ - ------------------i |
                                  \       225179981368524800   281474976710656000  /

        / 8575 \ /                / 1372 #2 \ \ /      27843348003288111           30791968342742367      \
     cos| ---- | | cosh(#2) - cosh| ------- | | | -------------------------- - -------------------------i |
        \  27  / \                \  15 pi  / / \ 36893488147419103232000000   2305843009213693952000000  /
   + ------------------------------------------------------------------------------------------------------
                                   /   2     769502681180271      20504211473148831   \
                          cosh(#2) | pi  + ------------------- - -------------------i |
                                   \       5629499534213120000   7036874417766400000  /

        / 12005 \ /                / 1372 #1 \ \ /        568231591903839            30791968342742367      \
     cos| ----- | | cosh(#1) - cosh| ------- | | | - ------------------------- + -------------------------i |
        \   27  / \                \  15 pi  / / \   2066035336255469780992000   6327233217282376204288000  /
   + --------------------------------------------------------------------------------------------------------
                                   /   2      769502681180271       20504211473148831   \
                          cosh(#1) | pi  + -------------------- - --------------------i |
                                   \       11033819087057715200   13792273858822144000  /

where

             /         2                                           \
             | 30625 pi      237500827524775     1265692066243755  |
   #1 == sqrt| --------- + ------------------ - -----------------i |
             \    1296     144115188075855872   36028797018963968  /

             /         2                                           \
             | 15625 pi      237500827524775     1265692066243755  |
   #2 == sqrt| --------- + ------------------ - -----------------i |
             \    1296     144115188075855872   36028797018963968  /

             /       2                                           \
             | 625 pi      237500827524775     1265692066243755  |
   #3 == sqrt| ------- + ------------------ - -----------------i |
             \   1296    144115188075855872   36028797018963968  /

             /       2                                           \
             | 625 pi      237500827524775     1265692066243755  |
   #4 == sqrt| ------- + ------------------ - -----------------i |
             \   144     144115188075855872   36028797018963968  /

             /       2                                           \
             | 625 pi      237500827524775     1265692066243755  |
   #5 == sqrt| ------- + ------------------ - -----------------i |
             \    16     144115188075855872   36028797018963968  /

[Dyuman Joshi]

Using pretty is not recommended.

If you want a numeric value as output, use a numeric data-type. And as you are dealing with large floating point numbers, double() would be the better choice.

"How can i define the calculation with infinity upper limit and odd number"

Let K=1 to Inf, Then 2*K-1 will be 1, 3, 5, .... Inf

Why is effective speed sound a complex value?

%% Analytical calculation

%% Input data
f = 15; % frequency [Hz]
omega = 2*pi.*f; % angular frequency
ho = 0.3*1e-3/2; % air gap height [m]
mu = 18.1e-6; % dynamic viscosity [Ns/m2]
rhoo = 1.206; % mean density [Kg/m3]
Co = 343; % undisturbed speed of sound [m/s]
Cp = 1004; % specific heat at constant pressure [J/KgK]
Cv = 716; % specific heat at constant volume [J/KgK]
lambda = 25.6e-3; % thermal conductivity [W/mK]
lx = 250*1e-3/2; % half length of plate [m]
ly = 180*1e-3/2; % half width of plate [m]
h = 0.05.*1e-3; % displacement amplitude
po = 1.01e5; % ambient pressure [N/m2]

%% Dimensionless Parameters
s = ho.*sqrt(rhoo.*omega./mu); % shear wave number
k = omega.*ho./Co; % reduced frequency
sigma = sqrt(mu.*Cp./lambda); % square root of Prandtl number
g = ho./lx; % nawroness of the gap
gamma = Cp/Cv; % ratio of specific heats
a = ly/lx; % aspect ratio of the panel
kx = omega*lx/Co; % normalised wave number along x-axis
ky = omega*ly/Co; % normalised wave number along y-axis
Bs =(tanh(s.*sqrt(i))./(s.*sqrt(i)))-1;
Bssigma =(tanh(s*sigma.*sqrt(i))./(s*sigma.*sqrt(i)))-1;
nssigma = (1+((gamma-1)/gamma).*Bssigma).^(-1);
Gamma = sqrt(gamma/(nssigma*Bs));
Ceffs = Co.*sqrt(nssigma.*Bs/gamma) % effective speed sound [m/s]

Ceffs = 43.3677 -45.4499i

syms x y q 

D = sqrt(((q*pi/(2*a))^2) - (kx*Gamma)^2);

SUMp=0;
a=1:2:9;
for jj=1:length(a)
    
     p = (2*nssigma*Gamma^2*kx*ky*h/pi).*symsum((((-1)^((q-1)/2)/(q*(D)^2)))*((cosh(D*x/kx)/cosh(D))-1)*cos(q*pi*y/(2*ky)),q,a(jj),a(jj));     

    SUMp=SUMp+p;
end

[Loïc M.]

Thanks, "double" is the correct function.

"Ceffs" doesn't apply in this equation.

Please could you provide me the correct line code to run with infinity upper limit and odd number?

Answer 1

[Walter Roberson]

Pi = sym(pi);
syms D x y kh kx ky

syms q Q
q = 2*Q - 1;
expression = (((-1)^((q-1)/2)/(q*(D)^2)))*((cosh(D*x/kx)/cosh(D))-1)*cos(q*Pi*y/(2*ky));
symsum(expression,Q,1,5)

ans = 

You might notice that the result is the same as if you had created the individual terms for odd q, and added the terms together. Which you could also do by taking the expression and subs() the q values in and sum() the result -- which is likely to be more efficient than forcing symsum() to analyze to see if it can find an expression and then ending up just adding the terms.

Can symsum() happen to find the infinite sum in this case? Let us test:

symsum(expression,Q,1,inf)

ans = 

Looks like Yes, at least under certain conditions.

[Loïc M.]

Thanks!

I'm trying now to plot the calculate function p(x,y) but i have an error.

%% Input data

f = 15; % frequency [Hz]

omega = 2*pi.*f; % angular frequency

ho = 0.3*1e-3/2; % air gap height [m]

mu = 18.1e-6; % dynamic viscosity [Ns/m2]

rhoo = 1.206; % mean density [Kg/m3]

Co = 343; % undisturbed speed of sound [m/s]

Cp = 1004; % specific heat at constant pressure [J/KgK]

Cv = 716; % specific heat at constant volume [J/KgK]

lambda = 25.6e-3; % thermal conductivity [W/mK]

lx = 250*1e-3/2; % half length of plate [m]

ly = 180*1e-3/2; % half width of plate [m]

h = ho*0.3; % displacement amplitude

po = 1.01e5; % ambient pressure [N/m2]

%% Dimensionless Parameters

s = ho.*sqrt(rhoo.*omega./mu); % shear wave number

k = omega.*ho./Co; % reduced frequency

sigma = sqrt(mu.*Cp./lambda); % square root of Prandtl number

g = ho./lx; % nawroness of the gap

gamma = Cp/Cv; % ratio of specific heats

a = ly/lx; % aspect ratio of the panel

kx = omega*lx/Co; % normalised wave number along x-axis

ky = omega*ly/Co; % normalised wave number along y-axis


Bs =(tanh(s.*sqrt(i))./(s.*sqrt(i)))-1;

Bssigma =(tanh(s*sigma.*sqrt(i))./(s*sigma.*sqrt(i)))-1;

nssigma = (1+((gamma-1)/gamma).*Bssigma).^(-1);

Gamma = sqrt(gamma/(nssigma*Bs));

syms x y q Q

D = sqrt(((q*pi/(2*a))^2) - (kx*Gamma)^2);

q = 2*Q-1;
     
p = (2*nssigma*Gamma^2*kx*ky*h/pi).*symsum((((-1)^((q-1)/2)/(q*(D)^2)))*((cosh(D*x/kx)/cosh(D))-1)*cos(q*pi*y/(2*ky)),Q,1,inf);     

P = double(p);

Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for variables.

Error in sym/double (line 729)
Xstr = mupadmex('symobj::double', S.s, 0);

x = linspace(-lx,lx);
y = linspace(-ly,ly);
[x,y] = meshgrid(x,y);
contourf(x,y,P);

[Walter Roberson]

It looks like you do not happen to meet the case where the summation has a closed form.

%% Input data

f = 15; % frequency [Hz]

omega = 2*pi.*f; % angular frequency

ho = 0.3*1e-3/2; % air gap height [m]

mu = 18.1e-6; % dynamic viscosity [Ns/m2]

rhoo = 1.206; % mean density [Kg/m3]

Co = 343; % undisturbed speed of sound [m/s]

Cp = 1004; % specific heat at constant pressure [J/KgK]

Cv = 716; % specific heat at constant volume [J/KgK]

lambda = 25.6e-3; % thermal conductivity [W/mK]

lx = 250*1e-3/2; % half length of plate [m]

ly = 180*1e-3/2; % half width of plate [m]

h = ho*0.3; % displacement amplitude

po = 1.01e5; % ambient pressure [N/m2]

%% Dimensionless Parameters

s = ho.*sqrt(rhoo.*omega./mu); % shear wave number

k = omega.*ho./Co; % reduced frequency

sigma = sqrt(mu.*Cp./lambda); % square root of Prandtl number

g = ho./lx; % nawroness of the gap

gamma = Cp/Cv; % ratio of specific heats

a = ly/lx; % aspect ratio of the panel

kx = omega*lx/Co; % normalised wave number along x-axis

ky = omega*ly/Co; % normalised wave number along y-axis


Bs =(tanh(s.*sqrt(i))./(s.*sqrt(i)))-1;

Bssigma =(tanh(s*sigma.*sqrt(i))./(s*sigma.*sqrt(i)))-1;

nssigma = (1+((gamma-1)/gamma).*Bssigma).^(-1);

Gamma = sqrt(gamma/(nssigma*Bs));

syms x y q Q

Pi = sym(pi);

q = 2*Q-1;
D = sqrt(((q*Pi/(2*a))^2) - (kx*Gamma)^2);

inner = (((-1)^((q-1)/2)/(q*(D)^2)))*((cosh(D*x/kx)/cosh(D))-1)*cos(q*Pi*y/(2*ky))

inner = 

p = (2*nssigma*Gamma^2*kx*ky*h/Pi).*symsum(inner,Q,1,20)

p = 

N = 5;
xvec = linspace(-lx,lx,N);
yvec = linspace(-ly,ly,N);
[X, Y] = meshgrid(xvec,yvec);

Ps = subs(p, {x,y}, {X, Y});

Ps(1,1)

ans = 

double(ans)

ans = -4.8260e+84 + 8.6220e+85i

If you use symsum to infinity then it cannot evalulate any results.

If you use symsum to 10 then the values are complex and on the order of 1e36. So this is diverging

%P = double(Ps);

%contourf(xvec, yec, P);

[Loïc M.]

Ok! How can i plot this similar contour with the previous equation of "p", please ?

I have this error

[Walter Roberson]

Your P includes the unresolved variable q

In your code version at https://www.mathworks.com/matlabcentral/answers/1934699-symsum-function-with-odd-number-and-infinity-upper-limit#comment_2675564 you had

syms x y q Q
D = sqrt(((q*pi/(2*a))^2) - (kx*Gamma)^2);
q = 2*Q-1;

Notice that uses q in D before you assign into q. Those two lines need to be reversed

syms x y q Q
q = 2*Q-1;
D = sqrt(((q*pi/(2*a))^2) - (kx*Gamma)^2);