Help for solving linear system that involves Bessel and Hankel Functions

I have tried starting with the following but I don't know how to proceed. Also, I don' t know to do it for the case of M different of N.

First I am defining the Bessel and Hankel Functions and their derivatives:

M= 10; 
N = M;
scaled = 1; % parameter for Bessel functions (to avoid overflow for large M)

for m = 0:M
    J_c(m+1) = besselj(m,k*r,scaled);
    H_s(m+1) = besselh(m,2,k*r,scaled);
    J_cp(m+1) = m*(besselj(m,k*r,scaled))./(k_c*a_c) - (besselj(m+1,k*r,scaled));     
    H_sp(m+1) = (1/2)*(besselh(m-1,2,k*r,scaled)-besselh(m+1,2,k*r,scaled)); 

Answer 1

[Torsten]

What kind of functions are the P_m^n and Q_m^n ?

Fix a finite upper bound N for the loops instead of Inf. Then you have a linear system of equations in A_0,...,A_N,B_0,...,B_N. Set up the coefficient matrix Mat and the right-hand side vector rhs and solve it using \

Use 2*N as the new upper bound for the loops instead of Inf and check whether the solutions of your first computation (with N as upper bound) don't differ much from those of this computation.

Example for solving a linear system in MATLAB:

Mat = [1 3 ; 4 6];
rhs = [-2;5];
sol = Mat\rhs

sol = 2×1

4.5000 -2.1667

[uranus]

The definition of Q and P.

[uranus]

I don't understand how exactly to solve the system with all As and all Bs. If it was a 2x2 it would be more straightforward

[Torsten]

If you write code to compute the coefficients of the linear system (14) and (15) using two matrices C1(N+1,N+1), C2(N+1,N+1) and two vectors of the right-hand sides B1(N+1), B2(N+1) (where N is the upper index for your two sums), I'll show you how to solve the system.

[uranus]

I have written the following code. I am not sure if it works well though.

( I have edited the code in the original comment)

clear all
close all
%N and M define the dimensions of the system
M = 10; 
N = 20;

scaled = 1; % parameter for Bessel functions (see matlab documentation)
eta = 1.5;  % frequency
phi = 0;

%Geometrical parameters  
a = 1 ; %  radius
multD = 1.5

multD = 1.5000

D = multD*a;%
beta=1;%m/s
omega = eta*pi*beta/a;
k=omega/beta;

H_p = zeros(M+1,1);
J_p = zeros(M+1,1);
B1 = zeros(M+1,1);
B2 = zeros(M+1,1);
C1 = zeros(M+1,N+1);
C2 = zeros(M+1,N+1);

for m = 0:M;
      J_p(m+1) = m*(besselj(m,k*a,scaled))./(k*a) - (besselj(m+1,k*a,scaled));     
      H_p(m+1) = (1/2)*(besselh(m-1,2,k*a,scaled)-besselh(m+1,2,k*a,scaled));       
      B1(m+1) = -2*(-1i^m+exp(-1i*2*k*D*cos(phi))*1i^m)*cos(m*phi);
      B2(m+1) = -2*(-1i^m-exp(-1i*2*k*D*cos(phi))*1i^m)*sin(m*phi);
      eps_m = [1,2*ones(1,M)];

    for n = 0:N  
       P(m+1,n+1) = (besselh(n+m,2,2*k*D,scaled) + (-1)^m*besselh(n-m,2,2*k*D,scaled));
       Q(m+1,n+1) = (besselh(n+m,2,2*k*D,scaled) - (-1)^m*besselh(n-m,2,2*k*D,scaled));
     end
end


for m = 0:M
  for n = 0:N
    if m==n
        C1(m+1,n+1) = (2/eps_m(m+1))*(H_p(m+1)/J_p(m+1)) + P(m+1,n+1);
    else
        C1(m+1, n+1) = P(m+1, n+1);
    end
  end
end

for m = 0:M
    C2(m+1,1) = 0;
end   

for m=1:M
  for n=1:N
    if m==n       
        C2(m+1,n+1) = (2/eps_m(m+1))*(H_p(m+1)/J_p(m+1)) + Q(m+1,n+1);
    else        
        C2(m+1, n+1) = Q(m+1, n+1);
    end
  end
end

sizeC2=size(C2)

sizeC2 = 1×2

11 21

  
A = C1\B1

A = 
   0.4505 + 1.0924i
  -0.4447 - 1.5784i
   0.2820 - 2.5257i
   0.5523 - 0.6733i
  -0.6121 + 0.1672i
   0.0000 + 0.0000i
   0.0062 + 0.5057i
   0.0000 + 0.0000i
   0.0912 - 0.0756i
   0.0522 - 0.0389i
   0.0186 - 0.0143i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0000 + 0.0000i
   0.0366 - 0.0318i
  -0.0111 + 0.0095i

B = C2\B2;

Warning: Rank deficient, rank = 10, tol = 2.673179e-08.

[Torsten]

Your first system reads

c00 * A0 + (a00 * A0 + a01 * A1 + ... + a0N * AN) = r0
c10 * A1 + (a10 * A0 + a11 * A1 + ... + a1N * AN) = r1
    ...
cN0 * AN + (aN0 * A0 + aN1 * A1 + ... + aNN * AN) = rN

your second system reads

d00 * B0 + (b00 * B0 + b01 * B1 + ... + b0N * BN) = s0
d10 * B1 + (b10 * B0 + b11 * B1 + ... + b1N * BN) = s1
...
dN0 * BN + (bN0 * B0 + bN1 * B1 + ... + bNN * BN) = sN

If you write both of them as one system, you get as ( 2*(N+1) x 2*(N+1) ) system matrix M

[(c00+a00)  a01 ...      a0N        0         0         ...      0 
  a10      (c10+a11) ... a1N        0         0         ...      0
               ......................................           
  aN0      aN1  ...      (cN0+aNN)  0         0         ...      0                 
  0        0    ...      0          (d00+b00) b01       ...      b0N
  0        0    ...      0          b10       (d10+b11) ...      b1N
               ......................................
  0        0    ...      0          bN0       bN1       ...     (dN0+bNN)]

and as (2*(N+1)x1) vector of the right-hand side

  rhs = [r0; r1; ... ;rN; s0; s1; ...; sN]

and you can solve for the (2*(N+1)x1) solution vector

AB = [A0 ;A1; ... ;AN; B0 ;B1; ... ;BN]

by

AB = M\rhs

The only thing left is to define the matrix coefficients and right-hand sides.

[uranus]

the for loops I have in my answer define the matrices C1, C2 and the right-hand sides B1, B2. Then I used the \ operator to solve the systems, as you decribe here. Doesn't that make my solution correct at least in terms of methodology?

[Torsten]

Doesn't that make my solution correct at least in terms of methodology?

No. The matrices must be square, and I don't know why you use two different dimensions m and n.

I suggest you first work without any loops, set N = 3, e.g., and fill in manually what the coefficients r00, ..., r30, d00,..., d30, a00,...,a33,b00,...,b33,r0,...,r3,s0,...,s3 are in this case. Then you can solve for A0,...,A3,B0,...,B3 as shown above.

[uranus]

In fact I did that, I wrote the first terms down and then constructed the matrices. It is from there that I saw that the matrices are not square and this is why I use two different dimensions n, m. So, I should check again my calculations I guess.

[Torsten]

If you determine the matrices and right-hand sides as shown above, you will end up with the correct problem size - either (N+1) x (N+1) for both A and B if you want to solve for A and B separately or 2*(N+1) x 2*(N+1) if you want to solve for A and B in one go.