Help for solving linear system that involves Bessel and Hankel Functions

2024 1 3
1 个回答
更新时间：2024 1 18
8 次查看（30 天）

0 个投票

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

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

回答（1 个）

Torsten
Torsten 2024-1-3
编辑：Torsten 2024-1-3
Ran in:

0 个投票

Ran in:
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
uranus 2024-1-10
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
Torsten 2024-1-10
编辑：Torsten 2024-1-10
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
uranus 2024-1-16
编辑：uranus 2024-1-16
Ran in:
Ran in:
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
Torsten 2024-1-16
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
uranus 2024-1-18
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
Torsten 2024-1-18
编辑：Torsten 2024-1-18
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
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
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.

请先登录，再进行评论。

类别

在 帮助中心 和 File Exchange 中查找有关 Special Functions 的更多信息

标签

提问:

uranus
uranus
2024-1-3

编辑：

Torsten
Torsten
2024-1-18

Community Treasure Hunt

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

Start Hunting!

Translated by