Help for solving linear system that involves Bessel and Hankel Functions

Question

uranus 2024-1-3

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2066056-help-for-solving-linear-system-that-involves-bessel-and-hankel-functions

编辑： Torsten 2024-1-18

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)); 

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-1-3

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2066056-help-for-solving-linear-system-that-involves-bessel-and-hankel-functions#answer_1382646

编辑：Torsten 2024-1-3

在 MATLAB Online 中打开

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

9 个评论
显示 7更早的评论隐藏 7更早的评论

uranus 2024-1-16

编辑：uranus 2024-1-16

在 MATLAB Online 中打开

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 2024-1-16

在 MATLAB Online 中打开

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 2024-1-18

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 2024-1-18

编辑：Torsten 2024-1-18

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.

请先登录，再进行评论。