T =
I have this system of linear equations and i need to solve it for unknowns. how i can solve it in MATLAB? we have 8 unknowns and 8 equations.
Mostra commenti meno recenti
$$
\begin{aligned}
& b_{0}^{l}=a_{0}^{l} A_{00}^{l}+\sum_{n=a_{n}}^{\infty} a_{n}^{l} A_{0} n^{\ell}+Z_{0}^{l} \\
& b_{m}^{l}=a_{0}^{l} A_{m 0}^{l}+\sum_{n=1}^{\infty} a_{n}^{l} A_{m n}^{l}+Z_{m}^{l} \\
& a_{0}^{l}=b_{0}^{l} B_{00}^{l}+\sum_{m=1}^{\infty} b_{m}{ }^{l} B_{0 m}^{l}+c_{0}^{l} C_{00}^{l}+\sum_{m=1}^{\infty} c_{m}^{l} C_{0 m}^{l}+Y_{0}^{l} \\
& a_{n}^{l}=b_{0}^{l} B_{n 0}^{l}+\sum_{m=1}^{\infty} b_{m}^{l} B_{n m}^{l}+c_{0}^{l} C_{n 0}^{l}+\sum_{m=1}^{\infty} c_{m}^{l} C_{n m}^{l}+Y_{n}^{l}
\end{aligned}
$$
and:
$$
\begin{aligned}
& c_{0}^{l}=d_{0}^{l} D_{00}^{\ell}+\sum_{n=1}^{\infty} d_{n}^{l} D_{0 n}^{l}+X_{0}^{l} \\
& c_{m}^{l}=d_{0}^{l} D_{m 0}^{\ell}+\sum_{n=1}^{\infty} d_{n}^{l} D_{m n}^{\ell}+X_{m}^{\ell} \\
& d_{0}^{\ell} = b_{0}^{l} E_{00}^{\ell}+\sum_{m=1}^{\infty} b_{m}^{\ell} E_{0 m}^{l}+c_{0}^{l} F_{00}^{l}+\sum_{m=1}^{\infty} c_{m}^{l} F_{0 m}^{l}+W_{0}^{l} \\
& d_{n}^{l}=b_{0}^{l} E_{n 0}^{l}+\sum_{m=1}^{\infty} b_{m}^{l} E_{n m}^{l}+c_{0}^{l} F_{n 0}^{l}+\sum_{m=1}^{\infty} c_{m}^{l} F_{n m}^{l}+W_{n}^{l}
\end{aligned}
$$
Unknowns are:-
$$
\begin{gathered}
\left(a_{0}^{l}, a_{n}^{l}\right) ;\left(b_{0}^{l}, b_{n}^{l}\right),\left(c_{0}^{l}, c_{n}^{l}\right) and \\
\left(d_{0}^{l}, d_{n}^{l}\right) .
\end{gathered}
$$
The above system of equations which is linear one is the latex code for these equations and i need to solve it for these unknowns. I am confused that how i can solved this system. how i can write these equations in matrix form 
to solved it for the unknowns? Help me.
to solved it for the unknowns? Help me.7 Commenti
Sam Chak
il 25 Mar 2025
Hi Javeria,
Conventionally, we can use the equationsToMatrix() function to transform a system of linear equations to matrix form, A·x = b. But I see many sums of infinite series.
javeria
il 25 Mar 2025
javeria
il 25 Mar 2025
Hi @javeria
If the series are convergent, you should see that the sum of its terms approaches a finite value as the number of terms increases. Have you evaluated all the series using the symsum() command? For example...
syms k x
T = ((-1)^k)*4/(2*k + 1)
F = symsum(T, k, 0, Inf)
javeria
il 25 Mar 2025
javeria
il 25 Mar 2025
Risposta accettata
javeria
il 14 Apr 2025
Più risposte (4)
Sam Chak
il 25 Mar 2025
0 voti
I see, @javeria. So, you would like to learn how to solve a system of linear equations. You can refer to the example at the following link and begin writing the code accordingly.
If you are able to create the equations in LaTeX, I believe you can certainly write those equations in MATLAB. You can also search for more examples on Google: "how to solve a system of linear equations with MATLAB".
1 Commento
javeria
il 25 Mar 2025
If you want to avoid writing your system as A*X = B, you can use a nonlinear solver like "fsolve". As the system is linear, it should converge in one step.
Otherwise - as @Sam Chak already said - you can use "equationsToMatrix" to let MATLAB compute A and B from your equations:
28 Commenti
@Torsten as you can see from the screen shoot these coefficients are defined like this and the rest also.
Then you will have to compute them first before applying "fsolve" or "equationsToMatrix":
M = 20;
N = 50;
for i = 1:M
for j = 1:N
A(i,j) = ...;
end
end
for j = 1:N
for i = 1:M
B(j,i) = ...;
end
end
And you don't have 8 equations for 8 unknowns. You have 2*(M+1) + 2*(N+1) equations for the unknowns a_0,a_1,...,a_N,b_0,b_1,...,b_M,c_0,c_1,...,c_M,d_0,d_1,...,d_N.
javeria
il 25 Mar 2025
In MATLAB, you must save the parameters in arrays with array indices increased by 1:
M = 20;
N = 50;
for i = 1:M+1
A(i,1) = ...; % this is A_(i-1),0 1 <= i <= M+1
end
for j = 2:N+1
A(1,j) = ...; % this is A_0,(j-1) 1 <= j <= N+1
end
for i = 2:M+1
for j = 2:N+1
A(i,j) = ...; % this is A_(i-1),(j-1) 2 <= i <= M+1, 2 <= j <= N+1
end
end
Sam Chak
il 26 Mar 2025
@javeria, I admit that at first, I did not examine your problem carefully due to the numerous subscripts and superscripts in the equations. My mind visually simplified what my eyes saw as the variables to be determined
and the rest are either coefficients in 
or constants in 
.
or constants in .In fact, this series alone 
already has infinitely many unknowns 
and n approaches ∞.
already has infinitely many unknowns and n approaches ∞.
javeria
il 26 Mar 2025
javeria
il 26 Mar 2025
@javeria, I am not very familiar with sequences. Otherwise, I'd have recognized the Bessel function. However, sharing your initial code and the mathematical context allows others to understand the specific problem you are facing. You will likely receive guidance on how to proceed, but generally, you will still need to write the code step-by-step as directed.
Click on the indentation icon 
to start defining the constants and sequences, or directly copy/paste from your MATLAB Editor, if you already have the code.
to start defining the constants and sequences, or directly copy/paste from your MATLAB Editor, if you already have the code.Also, avoid sharing information incrementally, as some individuals may lose motivation in the discussion if the initially suggested solution does not address the underlying problem adequately.
Hi @javeria
Upon computing the coefficients and constant terms, the primary challenge lies in solving a very large system of n linear equations. From your image, ideally, the series should represent the sum of infinitely many terms. However, you must truncate the series by selecting an appropriate value of n based on the desired magnitude of the remainder.
I conducted some simulations based on the number of variables in the linear system, and here are the results:
- Elapsed time is 7.818219 seconds for 100 variables.
- Elapsed time is 77.054258 seconds for 200 variables.
- Elapsed time is 352.742283 seconds for 300 variables.
- Elapsed time is 1104.857874 seconds for 400 variables.
- Elapsed time is 2743.320971 seconds for 500 variables.
tic
%% define the number of equations and variables
n = 100;
%% create symbolic variables x(1) to x(n)
syms x [1 n];
%% initialize an array to hold the equations
eqns= sym(zeros(1, n));
%% define the range for random integers
minValue = 1; % Minimum value
maxValue = 100; % Maximum value
%% generate random coefficients and constants
for i = 1:n
coeffs = randi([minValue, maxValue], 1, n); % random coefficients
constant= randi(n); % random constant term
eqns(i) = sum(coeffs.*x) == constant; % create the equation
end
%% convert to matrix form and solve the linear system
[A, b] = equationsToMatrix(eqns, x)
solx = vpa(linsolve(A, b))
toc
javeria
il 26 Mar 2025
@Sam Chak My problem involves many coefficients and values that must be predefined before solving the system. I am currently working on defining these values, which is where I am facing some difficulty. Once these coefficients and parameters are correctly defined, I will be able to proceed with solving the system step-by-step. I am looking for guidance on how to properly define these values and coefficients so that I can move forward with my solution. i need to define these coefficient matrix before solving the system. here this mode of bessel function i.e l is zero. The k_0 is find from dispersion equation for progressive mode and k_n from evenscent mode , the other parameter are defined as: \begin{itemize}
\item \( R_E = 1.0 \): The external radius of the cylinder.
\item \( h = 4 \cdot R_E \): The water depth, which is 4 times the external radius \( R_E \).
\item \( b = h - d = 2 \cdot R_E \): The distance from the bottom of the cylinder to the sea bed.
\item \( g = 9.81 \, \text{m/s}^2 \): The acceleration due to gravity.
\item \( N = 20 \): The number of terms used to truncate the series in regions E and I (external and internal regions).
\item \( M = 50 \): The number of terms used to truncate the series in region B (the region between the cylinder and the sea bed).
\item \( \text{ratios} = [1.031, 1.5, 2.0, 4.0] \): The ratios \( R_E / R_I \) that are to be tested, where \( R_I \) is the internal radius of the cylinder.
\item 
\end{itemize} Additionally, we consider the dispersion relations for both progressive and evanescent modes: \begin{itemize}
\item **Dispersion equation for progressive mode**:
\begin{equation*}
\omega^2 = k g tanh(kh)
\end{equation*}
\item **Dispersion equation for evanescent mode**: This describes the behavior of the wave in which the energy decays exponentially with distance from the boundary (i.e., these waves do not propagate but decay).
\begin{equation*}
\frac{\omega^2}{g} + k_n tan(k_n h) = 0
\end{equation*}
\end{itemize}
Torsten
il 26 Mar 2025
I am looking for guidance on how to properly define these values and coefficients so that I can move forward with my solution.
You expect guidance in determining physical constants ? In a MATLAB forum ?
javeria
il 27 Mar 2025
H_l(k_0*R_I) is besselh(l,k(1)*RI), K_l(k_n*RI) is besselk(l,k(n+1)*RI) ?
And here are the recursion formulas for the derivatives:
Why not doing it as I suggested: in three separate loops ?
A = zeros(M,N);
A(1,1) = ...; % Set A_00 according to formula 1
for j = 2:N
A(1,j) ...; % Set A_0n according to formula 2
end
for i = 2:M
A(i,1) = ...; % Set A_m0 according to formula 3
end
for i = 2:M
for j = 2:N
A(i,j) = ...; % Set A_mn according to formula 4
end
end
@Torsten @Sam Chak Okay thanks. i will define all these coeffiients before solving the system. So when i came to the linear equations solving part then I am little bit confuse that how i can now formulate my system of equations in matrix form that it take the previosly defined coefficients matrix correctly.
Sam Chak
il 28 Mar 2025
No @javeria, you do not need to define all these coefficients at once in your next post. Instead, make your post and code progressive. Focus on getting at least one part of the code working first for easy troubleshooting.
For example, define the coefficients for the 
series 
. Use small values for m and n to verify whether the suggested for-loop approach functions appropriately. If successful, then implement the code for the 
, 
, and 
series.
series . Use small values for m and n to verify whether the suggested for-loop approach functions appropriately. If successful, then implement the code for the , , and series.After you defined the matrix A and the vector Z,
Z = ...;
a = sym('a',[N 1]);
b = sym('b',[M 1]);
eqns12 = b == A*a + Z;
would define the first two lines of equations you included with a given M x 1 vector Z and unknown vectors a of size N x 1 and b of size M x 1 .
javeria
il 3 Apr 2025
% Define dimensions
M = 3; % For m = 1, 2, 3
N = 2; % For n = 1, 2
a = sym('a', [N 1]); % a1, a2
b = sym('b', [M 1]); % b1, b2, b3
c = sym('c', [M 1]); % c1, c2, c3
d = sym('d', [N 1]); % d1, d2
% Initialize coefficient matrices
A = zeros(M, N); % 4x3 (Aij)
B = zeros(N, M); % 3x4 (Bij)
C = zeros(N, M); % 3x4 (Cij)
D = zeros(M, N); % 4x3 (Dij)
E = zeros(N, M); % 3x4 (Eij)
F = zeros(N, M); % 3x4 (Fij)
% Define A (Aij)
A(1,1) = 2;
for j = 2:N
A(1,j) = 3;
end
for i = 2:M
A(i,1) = 4;
end
for i = 2:M
for j = 2:N
A(i,j) = 5;
end
end
% Define B (Bij)
B(1,1) = 1;
for j = 2:M
B(1,j) = 2;
end
for i = 2:N
B(i,1) = 3;
end
for i = 2:N
for j = 2:M
B(i,j) = 1;
end
end
% Define C (Cij)
C(1,1) = 2;
for j = 2:M
C(1,j) = 1;
end
for i = 2:N
C(i,1) = 0;
end
for i = 2:N
for j = 2:M
C(i,j) = 2;
end
end
% Define D (Dij)
D(1,1) = 1;
for j = 2:N
D(1,j) = 2;
end
for i = 2:M
D(i,1) = 3;
end
for i = 2:M
for j = 2:N
D(i,j) = 1;
end
end
% Define E (Eij)
E(1,1) = 1;
for j = 2:M
E(1,j) = 2;
end
for i = 2:N
E(i,1) = 3;
end
for i = 2:N
for j = 2:M
E(i,j) = 1;
end
end
% Define F (Fij)
F(1,1) = 2;
for j = 2:M
F(1,j) = 1;
end
for i = 2:N
F(i,1) = 0;
end
for i = 2:N
for j = 2:M
F(i,j) = 2;
end
end
% Define right-hand side vectors
Z = zeros(M, 1); % Z0, Z1, Z2, Z3
Y = zeros(N, 1); % Y0, Y1, Y2
X = zeros(M, 1); % X0, X1, X2, X3
W = zeros(N, 1); % W0, W1, W2
% Define Z (start at 1 and increment by 1)
Z(1) = 1;
for i = 2:M
Z(i) = 2;
end
% Define Y (start at 5 and increment by 1)
Y(1) = 5;
for i = 2:N
Y(i) = 9;
end
% Define X (start at 8 and increment by 1)
X(1) = 8;
for i = 2:M
X(i) = 9;
end
% Define W (start at 12 and increment by 1)
W(1) = 12;
for i = 2:N
W(i) = 10;
end
eqn1 = A*a-b+Z==0;
eqn2 = B*b+C*c-a+Y==0;
eqn3 = D*d-c+X==0;
eqn4 = E*b+F*c+W-d==0;
eqns = [eqn1;eqn2;eqn3;eqn4];
vars = [a;b;c;d];
[AA,bb] = equationsToMatrix(eqns,vars);
sol = AA\bb;
asol = sol(1:N)
bsol = sol(N+1:N+M)
csol = sol(N+M+1:N+M+M)
dsol = sol(N+M+M+1:N+M+M+N)
When N and M get larger, you should use
sol = double(AA)\double(bb);
asol = sol(1:N)
bsol = sol(N+1:N+M)
csol = sol(N+M+1:N+M+M)
dsol = sol(N+M+M+1:N+M+M+N)
javeria
il 3 Apr 2025
javeria
il 3 Apr 2025
Torsten
il 3 Apr 2025
Did you copy my full code ?
I changed
b = sym('b', [1 M]); % b1, b2, b3
a = sym('a', [1 N]); % a1, a2
c = sym('c', [1 M]); % c1, c2, c3
d = sym('d', [1 N]); % d1, d2
to
a = sym('a', [N 1]); % a1, a2
b = sym('b', [M 1]); % b1, b2, b3
c = sym('c', [M 1]); % c1, c2, c3
d = sym('d', [N 1]); % d1, d2
Maybe this is the problem.
javeria
il 3 Apr 2025
% Define dimensions
M = 4; % For m = 0, 1, 2, 3
N = 3; % For n = 0, 1, 2
% Initialize coefficient matrices
A = zeros(M, N); % 4x3 (Aij)
B = zeros(N, M); % 3x4 (Bij)
C = zeros(N, M); % 3x4 (Cij)
D = zeros(M, N); % 4x3 (Dij)
E = zeros(N, M); % 3x4 (Eij)
F = zeros(N, M); % 3x4 (Fij)
% Define A (Aij)
A(1,1) = 2;
for j = 2:N
A(1,j) = 3;
end
for i = 2:M
A(i,1) = 4;
end
for i = 2:M
for j = 2:N
A(i,j) = 5;
end
end
% Define B (Bij)
B(1,1) = 1;
for j = 2:M
B(1,j) = 2;
end
for i = 2:N
B(i,1) = 3;
end
for i = 2:N
for j = 2:M
B(i,j) = 1;
end
end
% Define C (Cij)
C(1,1) = 2;
for j = 2:M
C(1,j) = 1;
end
for i = 2:N
C(i,1) = 0;
end
for i = 2:N
for j = 2:M
C(i,j) = 2;
end
end
% Define D (Dij)
D(1,1) = 1;
for j = 2:N
D(1,j) = 2;
end
for i = 2:M
D(i,1) = 3;
end
for i = 2:M
for j = 2:N
D(i,j) = 1;
end
end
% Define E (Eij)
E(1,1) = 1;
for j = 2:M
E(1,j) = 2;
end
for i = 2:N
E(i,1) = 3;
end
for i = 2:N
for j = 2:M
E(i,j) = 1;
end
end
% Define F (Fij)
F(1,1) = 2;
for j = 2:M
F(1,j) = 1;
end
for i = 2:N
F(i,1) = 0;
end
for i = 2:N
for j = 2:M
F(i,j) = 2;
end
end
% Define right-hand side vectors
Z = zeros(M, 1); % Z0, Z1, Z2, Z3
Y = zeros(N, 1); % Y0, Y1, Y2
X = zeros(M, 1); % X0, X1, X2, X3
W = zeros(N, 1); % W0, W1, W2
% Define Z (start at 1 and increment by 1)
Z(1) = 1;
for i = 2:M
Z(i) = 2;
end
% Define Y (start at 5 and increment by 1)
Y(1) = 5;
for i = 2:N
Y(i) = 9;
end
% Define X (start at 8 and increment by 1)
X(1) = 8;
for i = 2:M
X(i) = 9;
end
% Define W (start at 12 and increment by 1)
W(1) = 12;
for i = 2:N
W(i) = 10;
end
% Define identity and zero submatrices
I1 = zeros(M, M); % 4x4 identity
I2 = zeros(N, N); % 3x3 identity
I3 = zeros(M, M); % 4x4 identity
I4 = zeros(N, N); % 3x3 identity
for i = 1:M
I1(i,i) = 1;
I3(i,i) = 1;
end
for i = 1:N
I2(i,i) = 1;
I4(i,i) = 1;
end
O1 = zeros(M, M); % 4x4 zero
O2 = zeros(M, N); % 4x3 zero
O3 = zeros(N, M); % 3x4 zero
O4 = zeros(N, N); % 3x3 zero
O5 = zeros(M, M); % 4x4 zero
O6 = zeros(M, N); % 4x3 zero
O7 = zeros(N, N); % 3x3 zero
% Assemble the full matrix G
G = [I1, -A, O1, O2;
-B, I2, -C, O4;
O5, O6, I3, -D;
-E, O7, -F, I4];
% Assemble the vector b
y = [Z; Y; X; W];
b = sym('b', [1 ,M]); % b0, b1, b2, b3
a = sym('a', [1, N]); % a0, a1, a2
c = sym('c', [1 ,M]); % c0, c1, c2, c3
d = sym('d', [1, N]); % d0, d1, d2
x = [b, a, c, d]'; % Vector of unknowns
% Solve the system
x_solution = inv(G) * y;
% Display all components
% Display results
disp('Matrix M (14x14):'); disp(M);
disp('Vector b (14x1):'); disp(b);
disp('Solution vector x:'); disp(x_solution);
% Extract bm, an, cm, dn from x_solution
bm = x_solution(1:M); % b0, b1, b2, b3
an = x_solution(M+1:M+N); % a0, a1, a2
cm = x_solution(M+N+1:M+N+M); % c0, c1, c2, c3
dn = x_solution(M+N+M+1:end); % d0, d1, d2
disp('bm (b0, b1, b2, b3):'); disp(bm);
disp('an (a0, a1, a2):'); disp(an);
disp('cm (c0, c1, c2, c3):'); disp(cm);
disp('dn (d0, d1, d2):'); disp(dn);
when i run my this code online it works but on my matlab it give me error i donot why
20 Commenti
% Define dimensions
M = 3; % For m = 0, 1, 2, 3
N = 2; % For n = 0, 1, 2
% Initialize coefficient matrices
A = zeros(M, N); % 4x3 (Aij)
B = zeros(N, M); % 3x4 (Bij)
C = zeros(N, M); % 3x4 (Cij)
D = zeros(M, N); % 4x3 (Dij)
E = zeros(N, M); % 3x4 (Eij)
F = zeros(N, M); % 3x4 (Fij)
% Define A (Aij)
A(1,1) = 2;
for j = 2:N
A(1,j) = 3;
end
for i = 2:M
A(i,1) = 4;
end
for i = 2:M
for j = 2:N
A(i,j) = 5;
end
end
% Define B (Bij)
B(1,1) = 1;
for j = 2:M
B(1,j) = 2;
end
for i = 2:N
B(i,1) = 3;
end
for i = 2:N
for j = 2:M
B(i,j) = 1;
end
end
% Define C (Cij)
C(1,1) = 2;
for j = 2:M
C(1,j) = 1;
end
for i = 2:N
C(i,1) = 0;
end
for i = 2:N
for j = 2:M
C(i,j) = 2;
end
end
% Define D (Dij)
D(1,1) = 1;
for j = 2:N
D(1,j) = 2;
end
for i = 2:M
D(i,1) = 3;
end
for i = 2:M
for j = 2:N
D(i,j) = 1;
end
end
% Define E (Eij)
E(1,1) = 1;
for j = 2:M
E(1,j) = 2;
end
for i = 2:N
E(i,1) = 3;
end
for i = 2:N
for j = 2:M
E(i,j) = 1;
end
end
% Define F (Fij)
F(1,1) = 2;
for j = 2:M
F(1,j) = 1;
end
for i = 2:N
F(i,1) = 0;
end
for i = 2:N
for j = 2:M
F(i,j) = 2;
end
end
% Define right-hand side vectors
Z = zeros(M, 1); % Z0, Z1, Z2, Z3
Y = zeros(N, 1); % Y0, Y1, Y2
X = zeros(M, 1); % X0, X1, X2, X3
W = zeros(N, 1); % W0, W1, W2
% Define Z (start at 1 and increment by 1)
Z(1) = 1;
for i = 2:M
Z(i) = 2;
end
% Define Y (start at 5 and increment by 1)
Y(1) = 5;
for i = 2:N
Y(i) = 9;
end
% Define X (start at 8 and increment by 1)
X(1) = 8;
for i = 2:M
X(i) = 9;
end
% Define W (start at 12 and increment by 1)
W(1) = 12;
for i = 2:N
W(i) = 10;
end
G = [A,-eye(M),zeros(M),zeros(M,N);...
-eye(N),B,C,zeros(N);...
zeros(M,N),zeros(M),-eye(M),D;...
zeros(N),E,F,-eye(N)];
rhs = [-Z;-Y;-X;-W];
% Solve the system
x_solution = G\rhs;
asol = x_solution(1:N)
bsol = x_solution(N+1:N+M)
csol = x_solution(N+M+1:N+M+M)
dsol = x_solution(N+M+M+1:N+M+M+N)
javeria
il 4 Apr 2025
how to fix this error
javeria
il 7 Apr 2025
Sam Chak
il 7 Apr 2025
Hi @javeria
The Symbolic Math Toolbox is included in the basic version of MATLAB Online, which can be used freely for 20 hours per month. You can test @Torsten's code there, and if both M and N are relatively small integers, it should work as long as it does not continuously compute for over 15 minutes. Additionally, MATLAB Online will automatically sign out after 15 minutes of idle time. These are the only drawbacks.
Even if you are away from your computer, you can also conveniently connect to MATLAB Online from your iPhone, iPad, or Android phone when you install the MATLAB Mobile app.
javeria
il 13 Apr 2025
javeria
il 13 Apr 2025
javeria
il 13 Apr 2025
next i want to make use of these values for matrix A and then plot result against k_0RE .
clc;
close all;
clear all;
tic;
N = 20; % Number of terms in region E and I
M= 50; % Number of term in region B
RE = 1.0; % Radius of external cylinder
ratios = [1.031,1.5, 2.0, 4.0]; % Ratios RE/RI you want to test
b = 2.0*RE; % The distance from seabed to the bottom of cylinder (h-d)=b
h = 4.0*RE; % Water depth.
g = 9.81; % Gravity acceleration
for r = 1:length(ratios)
RI = ratios(r) * RE; % Update RI based on the current ratio
%%%%%%%%%%%%%%%%%%%%%%%%%
%% finding the root \lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;
end
%%%%%%%%%%%%%%%%%%%%% Define matrix A %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A = zeros(M,N);
A(1,1) = cosh(k0*h)*sinh(k0*b)/(k0*b*(2*k0*h+sinh(2*k0*h)))*besselh(0., 1., k0*RE) / (k0 * (besselh(-1., 1., k0*RE))); % Set A_00
for j = 2:N
A(1,j) =cos(k(j)*h)*sin(k(j)*b) / (k(j)*b*(2*k(j)*h + sin(2*k(j)*h)))* besselk(0, k(j)*RE) / (-k(j)*besselk(-1, k(j)*RE)); % Set A_0n
end
for i = 2:M
A(i,1) =2 * cosh(k0*h) * (-1)^(i-1) * k0*sinh(k0*b) / (b * (2*k0*h + sinh(2*k0*h)) * (k0^2 + (m(i)^2))) * besselh(0, 1, k0*RE) / (k0*(besselh(-1, 1, k0*RE))); % Set A_m0
end
for i = 2:M
for j = 2:N
A(i,j) = 2 * (-1)^(i-1) * sin(k(j)*b) / (b * (2*k(j)*h + sin(2*k(j)*h)) * (k(j)^2 - m(i)^2)) * besselk(0, k(j)*RE) / (-k(j)*besselk(-1, k(j)*RE));% Set A_mn
end
end
end
like this just for refrence that i want my x axis like this
javeria
il 14 Apr 2025
javeria
il 14 Apr 2025
javeria
il 14 Apr 2025
If you change
for idx = 1:length(nu)
to
for idx = 1:1
and display e.g. the matrix A you get from
% Initialize matrix A
A = zeros(M, N);
% Set A_00 (first element)
A(1, 1) = cosh(k(1) * h) * sinh(k(1) * b) / (k(1) * b * (2 * k(1) * h + sinh(2 * k(1)* h))) * besselh(0, 1, k(1) * RE) / (k(1) * besselh(-1, 1, k(1) * RE));
% Set A_0n for j = 2:N
for j = 2:N
A(1, j) = cos(k(j) * h) * sin(k(j) * b) / (k(j) * b * (2 * k(j) * h + sin(2 * k(j) * h))) * besselk(0, k(j) * RE) / (-k(j) * besselk(-1, k(j) * RE));
end
% Set A_m0 for i = 2:M
for i = 2:M
A(i, 1) = 2 * cosh(k(1) * h) * (-1)^(i - 1) * k(1) * sinh(k(1) * b) / (b * (2 * k(1) * h + sinh(2 * k(1) * h)) * (k(1)^2 + m(i)^2)) * besselh(0, 1, k(1) * RE) / (k(1) * besselh(-1, 1, k(1)* RE));
end
% Set A_mn for i = 2:M, j = 2:N
for i = 2:M
for j = 2:N
A(i, j) = 2 * (-1)^(i - 1) * sin(k(j) * b) / (b * (2 * k(j) * h + sin(2 * k(j) * h)) * (k(j)^2 - m(i)^2)) * besselk(0, k(j) * RE) / (-k(j) * besselk(-1, k(j) * RE));
end
end
, you will see that almost all elements are NaN values. I didn't check further which term leads to the undefined values.
Further, you generate k as a 2-dimensional matrix, but you use k in the other equations as if it were a 1-dimensional vector k(j). This cannot be correct.
Did you pass the introductory course MATLAB Onramp to learn the basics of the language ? It is free of costs and will help you at least with the programming part:
javeria
il 15 Apr 2025
Hi @javeria
Mathematical operations such as 
and 
are known as indeterminate forms and can result in NaN (Not a Number). Furthermore, any arithmetic operation involving NaN will also yield NaN. By examining coef_A_0 from the bottom up, you can trace NaN all the way from c_1(1) to X_1. Since there is a pseudoinverse operation in X_1 involving the matrix H, it is essential to verify the computation of each element in H.
and are known as indeterminate forms and can result in NaN (Not a Number). Furthermore, any arithmetic operation involving NaN will also yield NaN. By examining coef_A_0 from the bottom up, you can trace NaN all the way from c_1(1) to X_1. Since there is a pseudoinverse operation in X_1 involving the matrix H, it is essential to verify the computation of each element in H.%% Testing javeria's code
tic;
% N = 20; % Number of terms in region E and I
% M = 50; % Number of term in region B
N = 2; % Number of terms in region E and I
M = 2; % Number of term in region B
RE = 1.0; % Radius of external cylinder
ratios = [1.031,1.5, 2.0, 4.0]; % Ratios RE/RI you want to test
b = 2.0*RE; % The distance from seabed to the bottom of cylinder (h-d)=b
h = 4.0*RE; % Water depth.
% nu = linspace(0.01, 10, 100);% Loop over the range of nu
nu = [0.01, 0.02];
%% finding the root \lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;
end
for r = 1:length(ratios)
RI = RE / ratios(r); % Update RI based on the current ratio
% coef_A_0_all = []; % To store results for different ratios
% Preallocate wave numbers k
k = zeros(length(nu), N); % Preallocate k matrix
% Compute wave numbers k using the algorithm of Chamberlin & Porter
for idx = 1:length(nu) % Loop over the range of nu
a2=nu(idx)*h;
p=a2*(1-4*(1-(1+a2)*exp(-2*a2))/(2*a2+sinh(2*a2)))^(-0.25); % approximation for k_0 (Eq. 19)
k(idx,1) = p; % Store k_0
% Compute evanescent roots k1, ..., kN-1 using Eq. 26 + Eq. 23
for i2=1:N-1
b2=i2*pi-pi/2*tanh(2*a2/(i2*pi*pi));
for j=1:2
u=b2/(1-2*(b2*tan(b2)+a2)*sin(2*b2)/(a2*(sin(2*b2)+2*b2)))^0.5;
b2=u;
end
k(idx, i2 + 1) = b2; % Store the result for k_n
end
% Convert all to dimensional values by dividing by h
k(idx, :) = k(idx, :) / h;
% Make k1...kN imaginary (evanescent)
k(idx, 2:N) = k(idx, 2:N) * complex(0, 1);
% Compute the matrix A for different ratios
% Initialize matrix A
A = zeros(M, N);
% Set A_00 (first element)
A(1, 1) = cosh(k(1) * h) * sinh(k(1) * b) / (k(1) * b * (2 * k(1) * h + sinh(2 * k(1)* h))) * besselh(0, 1, k(1) * RE) / (k(1) * besselh(-1, 1, k(1) * RE));
% Set A_0n for j = 2:N
for j = 2:N
A(1, j) = cos(k(j) * h) * sin(k(j) * b) / (k(j) * b * (2 * k(j) * h + sin(2 * k(j) * h))) * besselk(0, k(j) * RE) / (-k(j) * besselk(-1, k(j) * RE));
end
% Set A_m0 for i = 2:M
for i = 2:M
A(i, 1) = 2 * cosh(k(1) * h) * (-1)^(i - 1) * k(1) * sinh(k(1) * b) / (b * (2 * k(1) * h + sinh(2 * k(1) * h)) * (k(1)^2 + m(i)^2)) * besselh(0, 1, k(1) * RE) / (k(1) * besselh(-1, 1, k(1)* RE));
end
% Set A_mn for i = 2:M, j = 2:N
for i = 2:M
for j = 2:N
A(i, j) = 2 * (-1)^(i - 1) * sin(k(j) * b) / (b * (2 * k(j) * h + sin(2 * k(j) * h)) * (k(j)^2 - m(i)^2)) * besselk(0, k(j) * RE) / (-k(j) * besselk(-1, k(j) * RE));
end
end
%%%%%%%%%%%%%%%%%%%%%%%% DefineMatrix B %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
B=zeros(N,M);
B(1,1) = 4 * sinh(k(1)*b) / (RE * log(RE/RI) *cosh(k(1)*h)); %set B_00
for j=2:M
B(1,j) = m(j)*((besselk(0,m(j)*RI)*besseli(-1, m(j)*RE) + besseli(0, m(j)*RI)*besselk(-1,m(j)*RE)))...
/(besseli(0, m(j)*RE)*besselk(0,m(j)*RI) - besseli(0,m(j)*RI)*besselk(0,m(j)*RE))*4*k(1)^2*(-1)^(j-1)*sinh(k(1)*b)/(cosh(k(1)*h)*(k(1)^2+m(j)^2)); %Set B_0m
end
for i=2:N
B(i,1)=4*sin(k(i)*b)/(RE*log(RE/RI)*cosh(k(i)*h));% Set B_n0
end
for i=2:N
for j=2:M
B(i,j)= m(j)*((besselk(0,m(j)*RI)*besseli(-1, m(j)*RE) + besseli(0, m(j)*RI)*besselk(-1,m(j)*RE)))...
/(besseli(0, m(j)*RE)*besselk(0,m(j)*RI) - besseli(0,m(j)*RI)*besselk(0,m(j)*RE))*4*k(i)^2*(-1)^(j-1)*sin(k(i)*b)/(cos(k(i)*h)*(k(i)^2-m(j)^2)); % Set B_nm
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%Define Matrix C%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
C=zeros(N,M);
C(1,1) = (-4 * sinh(k(1) * b) )/ (RE*log(RE / RI) *cosh(k(1)*h)) ; %set C_00
for j=2:M
C(1,j)=-1/(besseli(0, m(j)*RE)*besselk(0,m(j)*RI)-besseli(0,m(j)*RI)*besselk(0,m(j)*RE))*4*k(1)^2*(-1)^(j-1)*sinh(k(1)*b)/(cosh(k(1)*h)*(k(1)^2+m(j)^2)); %Set C_0m
end
for i=2:N
C(i,1)= (-4 * sin(k(i) *b)) / (RE * log(RE/ RI)* cos(k(i) * h)); %Set C_n0
end
for i=2:N
for j=2:M
C(i,j)=-1/(RE*(besseli(0, m(j)*RE)*besselk(0,m(j)*RI)-besseli(0,m(j)*RI)*besselk(0,m(j)*RE)))*4*k(i)^2*(-1)^(j-1)*sin(k(i)*b)/(cos(k(i)*h)*(k(i)^2-m(j)^2)); % Set C_nm
end
end
%%%%%%%%%%%%%%%%Define matrix D %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
D = zeros(M, N);
D(1,1) = (cosh(k(1) * h) * cosh(k(1) * b)) / (k(1) * b * (2 * k(1) * h + sinh(2 * k(1) * h))) * (besselj(0, k(1) * RI)) / ((k(1)*besselj(-1, k(1) * RI)));%D_00
for j=2:N
D(1,j)= (cos(k(j) * h) * sin(k(j) * b) )/ (k(j) * b * (2 * k(j) * h + sin(2 * k(j) * h))) * (besseli(0, k(j) * RI) )/ ((k(j)*besseli(-1, k(j) * RI))); %D_0n
end
for i=2:M
D(i,1)=(2 *k(1)* (-1)^(i-1)* cosh(k(1) * h) * sinh(k(1) * b) )/ ( b * (2 * k(1) * h + sinh(2 * k(1)* h)) * (k(1)^2 + m(i)^2)) *(besselj(0, k(1) * RI)) / ((k(1)*besselj(-1, k(1) * RI))); % Set D_m0
end
for i=2:M
for j=2:N
D(i,j)= (2 *k(j)* (-1)^(i-1)* cos(k(j) * h) * sin(k(j) * b) )/ (b * (2 * k(j) * h + sin(2 * k(j)* h)) * (k(j)^2 - m(i)^2)) * (besseli(0, k(j) * RI)) / ((k(j)*besseli(-1, k(j) * RI)));% D_mn
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Define matrix E %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
E=zeros(N,M);
E(1,1)=4*sinh(k(1)*b)/(RI*log(RE/RI)*cosh(k(1)*h)); % E_00
for j=2:M
E(1,j)= 1./(RI*(besseli(0, m(j)*RE)*besselk(0,m(j)*RI) - besseli(0,m(j)*RI)*besselk(0,m(j)*RE)))*4*k(1)^2*(-1)^(j-1)*sinh(k(1)*b)/((k(1)^2+m(j)^2)*cosh(k(1)*h)); % E_0m
end
for i=2:N
E(i,1)=4*sin(k(i)*b)/(RI*log(RE/RI)*cos(k(i)*h)); % E_n0
end
for i=2:N
for j=2:M
E(i,j)= 1./(RI*(besseli(0, m(j)*RE)*besselk(0,m(j)*RI) - besseli(0,m(j)*RI)*besselk(0,m(j)*RE)))*4*k(i)^2*(-1)^(j-1)*sin(k(i)*b)/((k(i)^2-m(j)^2)*cos(k(i)*h)); % E_nm
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%Define matrix F %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%F%%%%%%%%%%%%%%%%%%%%%%
F=zeros(N,M);
F(1,1)=-4*sinh(k(1)*b)/(RI*log(RE/RI)*cosh(k(1)*h)); % F_00
for j=2:M
F(1,j)= -m(j)*(besseli(0,m(j)*RE)*besselk(-1, m(j)*RI) + besselk(0, m(j)*RE)*besseli(-1,m(j)*RI))...
/(besseli(0, m(j)*RE)*besselk(0,m(j)*RI) - besseli(0,m(j)*RI)*besselk(0,m(j)*RE))*4*k(1)^2*(-1)^(j-1)*sinh(k(1)*b)/((k(1)^2+m(j)^2)*cosh(k(1)*h)); %F_0m
end
for i=2:N
F(i,1)= -4*sin(k(i)*b)/(RI*log(RE/RI)*cos(k(i)*h)); % F_m0
end
for i=2:N
for j=2:M
F(i,j)=-m(j)*(besseli(0,m(j)*RE)*besselk(-1, m(j)*RI) + besselk(0, m(j)*RE)*besseli(-1,m(j)*RI))...
/(besseli(0, m(j)*RE)*besselk(0,m(j)*RI) - besseli(0,m(j)*RI)*besselk(0,m(j)*RE))*4*k(i)^2*(-1)^(j-1)*sin(k(i)*b)/((k(i)^2-m(j)^2)*cos(k(i)*h)); % F_nm
end
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % Define right-hand side vector G
G = zeros(2*M+2*N, 1); % Total length = M + N+ M+ N = 2M+2N
% Block 1: G (size M = 4)
for i = 1:M
if i == 1
G(i, 1) = (besselj(0, k(1)*RE) *sinh(k(1)*b))/((b*k(1))*cosh(k(1)*h)); % Z_0^l
else
G(i, 1) = (2*besselj(0, k(1)*RE) *k(1)*(-1)^(i-1)*sinh(k(1)*b))/(b*cosh(k(1)*h) *(k(1)^2 +m(i)^2)); %Z_m^l
end
end
% Block 2: Y (size N = 3)
for j = 1:N
if j == 1
G(j + M, 1) = -k(1)*besselj(-1, k(1)*RE) * (2*k(1)*h + sinh(2*k(1)*h)) /(cosh(k(1)*h)^2); % Y_0^l
else
G(j + M, 1) = 0; % Y_n^l
end
end
% Block 3: X (size M = 4)
for i = 1:M
G(i + M + N, 1) = 0; % X_0^l, X_m^l
end
% Block 4: W (size N = 3)
for j = 1:N
G(j + M + N + M, 1) = 0; % W_0^l, W_n^l
end
% Define identity and zero submatrices
I1 = zeros(M, M); % MxM identity
I2 = zeros(N, N); % NxN identity
I3 = zeros(M, M); % MxM identity
I4 = zeros(N, N); % N*N identity
for i = 1:M
I1(i,i) = 1;
I3(i,i) = 1;
end
for i = 1:N
I2(i,i) = 1;
I4(i,i) = 1;
end
%%%%%%%%%%%% Define zero submatrices%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
O1 = zeros(M, M); %
O2 = zeros(M, N); %
O3 = zeros(N, N); %
O4 = zeros(M, M); %
O5 = zeros(M, N); %
O6 = zeros(N, N); %
% Construct the matrix H while solving H X = G
H = [I1, -A, O1, O2; ...
-B, I2, -C, O3; ...
O4, O5, I3, -D; ...
-E, O6, -F, I4]
% % Solution of system of equations using inv or pinv (pscheudo inverse)
% we can use X_1=pinv(H)*G
X_1 = pinv(H)*G
% Extract the unknowns from the solution vector X_1
b_1 = zeros(1, M); % b0, b1, b2, b3
a_1 = zeros(1, N); % a0, a1, a2
c_1 = zeros(1, M); % c0, c1, c2, c3
d_1 = zeros(1, N); % d0, d1, d2
% Block 1: b_1 (size M = 4)
for i = 1:M
b_1(i) = X_1(i, 1);
end
% Block 2: a_1 (size N = 3)
for j = 1:N
a_1(j) = X_1(j + M, 1);
end
% Block 3: c_1 (size M = 4)
for i = 1:M
c_1(i) = X_1(i + M + N, 1);
end
% Block 4: d_1 (size N = 3)
for j = 1:N
d_1(j) = X_1(j + M + N + M, 1);
end
term1(idx) = c_1(1);
coef_A_0(idx) = abs( ( c_1(1)*cosh(k(1)*h)^2 )/( 2*k(1)*h + sinh(2*k(1)*h) ) * k(1)*besselj(-1, k(1)*RI) );
end
end
toc;
disp(term1) % a term "c_1(1)" in coef_A_0
javeria
il 15 Apr 2025
Torsten
il 15 Apr 2025
% % This is the code for finding roots of the dispersion relation
% % K = k * tanh(kh)
% % here K = nu^2/g, nu is angular frequency and h is water depth
% function [ k1 ] = dispersion( nu,h,N)
RE = 1.0; % Radius of external cylinder
ratios = [1.031, 1.5, 2.0, 4.0]; % Ratios RE/RI you want to test
b = 2.0 * RE; % Distance from seabed to the bottom of cylinder (h-d)=b
h = 4.0 * RE; % Water depth
N=3; % how many roots you need? Change your N to get roots
M = 4;% Change your M according to u
Nu = [0.1;0.2]; % angular frequency
h=4.0; % water depth
%% finding the root \lemda_m =m*pi/b%%%%%%
for j=1:M
m(j)=(pi*(j-1))/b;
end
for nn = 1:numel(Nu)
nu = Nu(nn);
k = zeros(1, N);
% % Use the algorithm of Chamberlin & Porter
a2=nu*h;
p=a2*(1-4*(1-(1+a2)*exp(-2*a2))/(2*a2+sinh(2*a2)))^(-0.25); % approximation for k_0 (Eq. 19)
k(1)=p;
for i2=1:N-1
b2=i2*pi-pi/2*tanh(2*a2/(i2*pi*pi));%aproximation for k_n (Eq. 26)
for j=1:2
u=b2/(1 - (2*(b2*tan(b2) + a2)*sin(2*b2))/(a2*(sin(2*b2)+2*b2)))^0.5; % From article (Eq. 23), an upper bound for k_n, iterated to refine the roots
b2=u;
end
k(i2+1)=b2;
end
k=k/h;
k(:,2:N)=k(:,2:N)*complex(0,1); % print all roots where first one is real and remaining ones are purely imaginary
for r = 1:length(ratios)
RI = RE / ratios(r); % Update RI based on the current ratio
% Example: Initialize matrix A as a zero matrix of size NxN
A = zeros(M, N);
% Set A_00 (first element)
A(1, 1) = cosh(k(1) * h) * sinh(k(1) * b) / (k(1) * b * (2*k(1)*h + sinh(2 * k(1)* h))) * besselh(0, 1, k(1) * RE) / (k(1) * besselh(-1, 1, k(1) * RE));
% Set A_0n for j = 2:N
for j = 2:N
A(1, j) = cos(k(j) * h) * sin(k(j) * b) / (k(j) * b * (2*k(j)*h + sin(2*k(j) * h))) * besselk(0, k(j) * RE) / (-k(j) * besselk(-1, k(j) * RE));
end
% Set A_m0 for i = 2:M
for i = 2:M
A(i, 1) = 2 * cosh(k(1) * h) * (-1)^(i - 1) * k(1) * sinh(k(1) * b) / (b*(2 * k(1) * h + sinh(2 * k(1) * h)) * (k(1)^2 + m(i)^2)) * besselh(0, 1, k(1) * RE) / (k(1) * besselh(-1, 1, k(1)* RE));
end
% Set A_mn for i = 2:M, j = 2:N
for i = 2:M
for j = 2:N
A(i, j) = 2 * (-1)^(i - 1) * sin(k(j) * b) / (b * (2 * k(j) * h + sin(2 * k(j) * h)) * (k(j)^2 - m(i)^2)) * besselk(0, k(j) * RE) / (-k(j) * besselk(-1, k(j) * RE));
end
end
% After calculating matrix A, you can process or print it as required
fprintf('Matrix A for ratio %.3f and RI %.3f:\n', ratios(r), RI);
disp(A);
end
end
javeria
il 15 Apr 2025
javeria
il 16 Apr 2025
0 voti
Categorie
Scopri di più su Mathematics in Centro assistenza e File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
