I want to chek my function for solving system of equaton

Hi

i write matlba function for solve a systme from a EM problem

function [Is] = current()
    [f,N,Nc,a1,a2,ra,k0,Z0,lambda ,gap_angle] = parameter();
    
    % Angular positions φ_n
    n=1:N;
    phi_n = -gap_angle+ ((2*n - 1) * pi / N);  % Compute φ_n
    
    % Initialize matrices
    A = zeros(4*N, 2*N);          % Coefficient matrix (4N equations, 2N unknowns)
    b = zeros(4*N, 1);            % Right-hand side vector
    
    % Fill the coefficient matrix A and vector b
    for m = 1:N
        for n = 1:N
            % Compute angular difference and distances R1 and R2
            phi_diff = phi_n(m) - phi_n(n);
            R1 = sqrt(ra^2 + a1^2 - 2*ra*a1*cos(phi_diff));
            R2 = sqrt(ra^2 + a2^2 - 2*ra*a2*cos(phi_diff));
            
            % Bessel functions
            H0_R1 = besselh(0, 2, k0*R1);  % H_0^(2)(k0*R1)
            H0_R2 = besselh(0, 2, k0*R2);  % H_0^(2)(k0*R2)
            H1_R1 = besselh(1, 2, k0*R1);  % H_1^(2)(k0*R1)
            H1_R2 = besselh(1, 2, k0*R2);  % H_1^(2)(k0*R2)
            
            % Equation (5): Rows 1:N
            A(m, n)=(k0*Z0/4) * H0_R1;
            
            % Equation (6): Rows N+1:2N
            A(m+N, n+N) = H0_R2;
            
            % Equation (7): Rows 2N+1:3N
            A(m+2*N, n)   = (k0*Z0/4) * H0_R1;
            A(m+2*N, n+N) = -(k0*Z0/4) * H0_R2;
            
            % Equation (8): Rows 3N+1:4N
            A(m+3*N, n)   = (1j*k0/4) * ((ra - a1*cos(phi_diff)) / R1) * H1_R1;
            A(m+3*N, n+N) = -(1j*k0/4) * ((ra - a2*cos(phi_diff)) / R2) * H1_R2;
        end
        
        % Right-hand side vector b
        b(m)       =E_inc_z(ra, phi_n(m));     % Incident E_z for Eq (5)
        b(m+N)     = 0;                        % Zero for Eq (6)
        b(m+2*N)   =E_inc_z(ra, phi_n(m));     % Incident E_z for Eq (7)
        b(m+3*N)   = H_inc_phi(ra, phi_n(m));   % Incident H_phi for Eq (8)
    end
    
    % Solve the system using linsolve
    x = linsolve(A, b);
    % Extract unknowns w_n^(1) and w_n^(2)
    w1 = x(1:N);
    w2 = x(N+1:2*N);
    Is = [w1; w2]
    % Display results
    %disp('w_n^(1):');
    %disp(w1);
    %disp('w_n^(2):');
    %disp(w2);

end    % added by Sam (Editor)

the system is in

the function is correct ?

thank

George

Answer 1

[Parag]

Hi, I saw you code and executed MATLAB code in R2024a and I find it correct.

function [Is] = current()

[f,N,Nc,a1,a2,ra,k0,Z0,lambda ,gap_angle] = parameter();

% Angular positions φ_n

n=1:N;

phi_n = -gap_angle+ ((2*n - 1) * pi / N); % Compute φ_n

% Initialize matrices

A = zeros(4*N, 2*N); % Coefficient matrix (4N equations, 2N unknowns)

b = zeros(4*N, 1); % Right-hand side vector

% Fill the coefficient matrix A and vector b

for m = 1:N

for n = 1:N

% Compute angular difference and distances R1 and R2

phi_diff = phi_n(m) - phi_n(n);

R1 = sqrt(ra^2 + a1^2 - 2*ra*a1*cos(phi_diff));

R2 = sqrt(ra^2 + a2^2 - 2*ra*a2*cos(phi_diff));

% Bessel functions

H0_R1 = besselh(0, 2, k0*R1); % H_0^(2)(k0*R1)

H0_R2 = besselh(0, 2, k0*R2); % H_0^(2)(k0*R2)

H1_R1 = besselh(1, 2, k0*R1); % H_1^(2)(k0*R1)

H1_R2 = besselh(1, 2, k0*R2); % H_1^(2)(k0*R2)

% Equation (5): Rows 1:N

A(m, n)=(k0*Z0/4) * H0_R1;

% Equation (6): Rows N+1:2N

A(m+N, n+N) = H0_R2;

% Equation (7): Rows 2N+1:3N

A(m+2*N, n) = (k0*Z0/4) * H0_R1;

A(m+2*N, n+N) = -(k0*Z0/4) * H0_R2;

% Equation (8): Rows 3N+1:4N

A(m+3*N, n) = (1j*k0/4) * ((ra - a1*cos(phi_diff)) / R1) * H1_R1;

A(m+3*N, n+N) = -(1j*k0/4) * ((ra - a2*cos(phi_diff)) / R2) * H1_R2;

end

% Right-hand side vector b

b(m) =E_inc_z(ra, phi_n(m)); % Incident E_z for Eq (5)

b(m+N) = 0; % Zero for Eq (6)

b(m+2*N) =E_inc_z(ra, phi_n(m)); % Incident E_z for Eq (7)

b(m+3*N) = H_inc_phi(ra, phi_n(m)); % Incident H_phi for Eq (8)

end

% Solve the system using linsolve

x = linsolve(A, b);

% Extract unknowns w_n^(1) and w_n^(2)

w1 = x(1:N);

w2 = x(N+1:2*N);

Is = [w1; w2]

% Display results

%disp('w_n^(1):');

%disp(w1);

%disp('w_n^(2):');

%disp(w2);

end % added by Sam (Editor)

function [f, N, Nc, a1, a2, ra, k0, Z0, lambda, gap_angle] = parameter()

% Define or load the parameters here

f = 10e9; % Frequency (Hz)

N = 10; % Number of elements (example)

Nc = 4; % Some parameter Nc (example)

a1 = 0.5; % Radius 1 (example)

a2 = 0.6; % Radius 2 (example)

ra = 0.8; % Radius a (example)

k0 = 2 * pi * f / 3e8; % Wave number (example with speed of light)

Z0 = 377; % Impedance of free space (Ohms)

lambda = 3e8 / f; % Wavelength (meters)

gap_angle = pi/4; % Gap angle (example)

end

function E = E_inc_z(ra, phi)

% Define the incident electric field E_z (for example, a plane wave)

E = exp(-1j * ra * cos(phi)); % Example expression, modify as per your requirements

end

function H = H_inc_phi(ra, phi)

% Define the incident magnetic field H_phi (for example, a plane wave)

H = exp(-1j * ra * sin(phi)); % Example expression, modify as per your requirements

end

Is = current();

disp(Is);

Below are the key points for further explanation:

• the result is a 20×1 matrix, which matches the output size expecting for w1 and w2 combined (since w1 and w2 are each of size N×1N and N= 10 in setup).
• The first 10 values (from 1 to 10) are likely associated with the first unknowns (w1), and the second 10 values (from 11 to 20) are likely associated with the second unknowns (w2).
• For equation 5, the term (k0*Z0/4) * H0_R1 is correctly placed in the matrix A for the first set of equations. The right-hand side b(m) is set to E_inc_z(ra, phi_n(m)), which matches the equation.
• For equation 6, The term H0_R2 is correctly placed in the matrix A for the second set of equations. The right-hand side b(m+N) is set to 0, which matches the equation.
• For equation 7, The terms (k0*Z0/4) * H0_R1 and -(k0*Z0/4) * H0_R2 are correctly placed in the matrix A. The right-hand side b(m+2*N) is set to E_inc_z(ra, phi_n(m)), which matches the equation.
• For equation 8, The terms (1j*k0/4) * ((ra - a1*cos(phi_diff)) / R1) * H1_R1 and -(1j*k0/4) * ((ra - a2*cos(phi_diff)) / R2) * H1_R2 are correctly placed in the matrix A. The right-hand side b(m+3*N) is set to H_inc_phi(ra, phi_n(m)), which matches the equation.

Answer 2

[george veropoulos]

a better version is

function [w1,w2] = currentMAS()

[~,N,a1,a2,ra,k0,Z0,~ ,gap_angle] = parameter();

% Angular positions φ_n

phi_n = -gap_angle + (2*(1:N) - 1) * ( pi/ N);

phi_m = phi_n; % Matching points are the same as source points

% Compute distances R1_nm and R2_nm

R1_nm = sqrt(ra^2 + a1^2 - 2 * ra * a1 * cos(phi_m' - phi_n));

R2_nm = sqrt(ra^2 + a2^2 - 2 * ra * a2 * cos(phi_m' - phi_n));

% Compute Hankel functions

H0_2_R1 = besselh(0, 2, k0 * R1_nm);

H0_2_R2 = besselh(0, 2, k0 * R2_nm);

H1_2_R1 = besselh(1, 2, k0 * R1_nm);

H1_2_R2 = besselh(1, 2, k0 * R2_nm);

% Preallocate coefficient matrix A and right-hand side B

A = zeros(2*N, 2*N); % Adjust size based on number of unknowns

B = zeros(2*N, 1);

% Construct the system based on conditions

for m = 1:N

if abs(phi_m(m)) < gap_angle

% Equations (5) and (6) for |phi_m| < gamma

A(m, 1:N) = (k0 *Z0 / 4) * H0_2_R1(m, :); % w_n^(1)

A(m, N+1:end) = 0; % w_n^(2) terms

B(m) = E_inc_z(ra, phi_m(m)); % Incident field

A(N+m, 1:N) = 0; % w_n^(1) terms

A(N+m, N+1:end) = H0_2_R2(m, :); % w_n^(2)

B(N+m) = 0;

elseif abs(phi_m(m)) > gap_angle

% Equations (7) and (8) for |phi_m| > gamma

A(m, 1:N) = (k0 * Z0 / 4) .* H0_2_R1(m, :); % w_n^(1)

A(m, N+1:end) = -(k0 * Z0 / 4).* H0_2_R2(m, :); % w_n^(2)

B(m) =E_inc_z(ra, phi_m(m));

A(N+m, 1:N) = (1j * k0 / 4) * H1_2_R1(m, :) .* (ra - a1 * cos(phi_m(m) - phi_n)) ./ R1_nm(m, :);

A(N+m, N+1:end) = -(1j * k0 / 4) * H1_2_R2(m, :) .* (ra - a2 * cos(phi_m(m) - phi_n)) ./ R2_nm(m, :);

B(N+m) =H_inc_phi(ra, phi_m(m));

end

end

% Solve the system

cond_A = cond(A);

disp(cond_A);

w = linsolve(A, B);

%w = A\B;

% Extract solutions

w1 = w(1:N); % w_n^(1)

w2 = w(N+1:end); % w_n^(2)

end