How to fix the code for the derivative of the Hankel function asymtotic expansion? [I added my code here]

I want to write the prime of this hankel function expansion, where , zeta, B0, ka are given, nu is complex which is also given in the question, How can I write the derivatiove of this Hankel function first kind equation 21? I am providing here the equations in attached file also I am providing the code I wrote, I am not sure this is the exact way to write the prime or not. Can anybody help me to fix this derivative part of this ciode?

%% Constants
mu0 = 4 * pi * 1e-7;       % Permeability of free space
eps0 = 8.854e-12;          % Permittivity of free space
c = 3e8;                   % Speed of light
freq = [1e9];                % Frequency (1 GHz)
lambda0 = c / (freq);        % Free-space wavelength
k0 = 2 * pi ./ (lambda0);     % Free-space wavenumber
b = 20 ./ (k0);               % Outer radius of the dielectric coating (~0.955 meters)
eps_r = 4;  
%mu_r= mu/mu0% Relative permittivity
eps1 = eps_r * eps0;       % Permittivity of dielectric
k1 = k0 * sqrt(eps_r);  
xl = -(2.33810741);
omega= 2*pi*freq;
d = 0.30*lambda0;% thickness of the dielectric material in cylider

a = b - d; % Compute inner radius of the cylinder

%% Compute k1*a and k1*b
k1a = k1 * (a);
k1b = k1 * (b);
% 
%disp (k1a)
% %k1= omega* sqrt (eps1); %% if Naishadham's paper 
%% 
%% nu_equation for k1a case
term1 = k1 * a;
term2 = - exp(1i * pi/3) * xl * (k1*a)^(1/3);
% Term 3 
term3 = (1/60) ...
    * exp(1i * (2*pi/3)) ...
    * (xl^2) ...
    * ((1/2)*k1*a)^(-1/3);
% Term 4:
term4 = - ((xl^3 + 10) / 1400) ...
    * ((1/2)*k1*a)^(-1);
% Term 5:
term5 = - exp(1i * pi/3) ...
    * (281*xl^4 + 10440*xl) / 4536000 ...
    * ((1/2)*k1*a)^(-5/3);
% Term 6:
term6 = - exp(1i * 2*pi/3) ...
    * (73769*xl^5 + 6624900*xl^2) / 10478160000 ...
    * ((1/2)*k1*a)^(-7/3);
% Term 7:
term7 = ( (93617*xl^6 + 16495400*xl^3 - 1744600) / 100900800000 ) ...
    * ((1/2)*k1*a)^(-3);
% Term 8 (Big-O term):
term8 = (k1*a)^(-11/3);
nu_k1a = term1 + term2 + term3+term4+term5+term6+term7+term8;
disp(nu_k1a);

  40.0812 + 6.7302i

%% write expression for zeta_k1a using equation 16
bracketTerm = ...
    log( ( nu_k1a + sqrt(nu_k1a^2 - k1a^2 ) ) / k1a ) ...
    - sqrt( (nu_k1a^2 - k1a^2) / (nu_k1a^2) );
zeta_k1a = ( (3/2)^(2/3) ) * ( bracketTerm )^(2/3);
disp('zeta_k1a =');

zeta_k1a =

disp(zeta_k1a);

   0.1480 + 0.2001i

%% write beta0 zeta_k1a case
B0_zeta_k1a = -5/(48 * zeta_k1a^2) ...
    + (1 / (24 * sqrt(zeta_k1a))) * ...
    ( (5 * nu_k1a^3) / ((nu_k1a^2 - (k1*a)^2)^(3/2)) ...
    - (3 * nu_k1a)   / ((nu_k1a^2 - (k1*a)^2)^(1/2)) );
disp('B0_zeta_k1a');

B0_zeta_k1a

disp(B0_zeta_k1a);

   0.0193 + 0.0019i

%% Hankel function first kind expansion k1a case 
% 2) Airy argument
% Now we have nu_k1a in place of nu, and zeta_k1a in place of zeta.
zArg = exp(1i * (2*pi/3)) * (nu_k1a^(2/3)) * zeta_k1a;
% 3) Compute Airy function and derivative
% airy(0,z) = Ai(z)
% airy(1,z) = Ai'(z)
AiVal      = airy(0, zArg);  % Ai(zArg)
AiPrimeVal = airy(1, zArg);  % Ai'(zArg)
% prefactor = 2*sqrt(2) * [ (nu_k1a^2 * zeta_k1a) / (nu_k1a^2 - (k1a)^2 ) ]^(1/4)
prefactor = 2*sqrt(2) * ...
    ( (nu_k1a^2 * zeta_k1a) / (nu_k1a^2 - (k1a)^2 ) )^(1/4);
% bracket = e^(-i*pi/3)*AiVal*nu_k1a^(-1/3)
%         + e^( i*pi/3)*AiPrimeVal*nu_k1a^(-5/3)*B0_zeta_k1a
term1 = exp(-1i * pi/3) ...
    * AiVal ...
    * (nu_k1a^(-1/3));
term2 = exp(1i * pi/3) ...
    * AiPrimeVal ...
    * (nu_k1a^(-5/3)) ...
    * B0_zeta_k1a;
bracket = term1 + term2;
Term3 = 1 + 1/(nu_k1a^2);
% ========================================
H_nu_k1a_approx = prefactor * bracket * Term3;
% ======================
% 8) Display the result
% ======================
disp('Asymptotic Hankel H_nu^(1)(k1*a) = ');

Asymptotic Hankel H_nu^(1)(k1*a) = 

disp(H_nu_k1a_approx);

  -0.1275 + 0.2372i

%% write the derivative of Hankel function first kind k1a case
%%%% [derivative of equation 21] 
%2) Define the Hankel function (first kind) as a function handle
%Hfun(a)
%2) Define the big-O replacement factor:
Term3 = 1 + 1/(nu_k1a^2);
% 3) Build a function handle for H_nu^(1)(k1*a)
%    We implement the prefactor, bracket, and Term3 as described above.
Hfun = @(a) ...
    ( 2*sqrt(2) ...
    * ( (nu_k1a^2 * zeta_k1a) ./ (nu_k1a^2 - (k1.*a).^2 ) ).^(1/4) ) ...
    .* ( ...
    exp(-1i*pi/3)*AiVal            * nu_k1a^(-1/3) ...
    + exp( 1i*pi/3)*AiPrimeVal       * nu_k1a^(-5/3) * B0_zeta_k1a ...
    ) ...
    .* Term3;
% Explanation of the main pieces:
%  - The prefactor = 2*sqrt(2)* [ (nu_k1a^2*zeta_k1a)/(nu_k1a^2 - (k1*a)^2 ) ]^(1/4).
%  - The bracket   = e^(-i*pi/3)*AiVal*nu_k1a^(-1/3)
%                    + e^( i*pi/3)*AiPrimeVal*nu_k1a^(-5/3)*Beta0_k1a.
%  - We multiply by (1 + 1/nu_k1a^2) for the big-O term.
% 4) Finite-difference derivative at a
da = 1e-6;  % small step
H_plus  = Hfun(a + da);
H_minus = Hfun(a - da);
H_prime_approx = (H_plus - H_minus) / (2*da);
disp('Approx derivative d/da [H_nu^(1)(k1*a)] at a :');

Approx derivative d/da [H_nu^(1)(k1*a)] at a :

disp(H_prime_approx)

   0.2070 + 0.2747i

[Torsten]

Do you have reason to believe your finite-difference approximation of the derivative is wrong ?

[SSWH2]

Yes, because when I am writing my derivative like this way mentioned bellow- I am getting different answer, so I don't know which approach of writing this equation derivative is correct..

%% Write the derivative of Hankel function first kind k1a case
% 1) Save the original a and k1a
a_orig   = a;

Unrecognized function or variable 'a'.

k1a_orig = k1a;
% 2) Choose a small step for finite difference
da = 1e-6;
% ==========================
% Evaluate at a_plus = a + da
% ==========================
a_plus   = a_orig + da;
k1a_plus = k1 * a_plus;
% --- Now replicate your lines for nu_k1a, zeta_k1a, B0_zeta_k1a, H_nu_k1a_approx ---
%     but with (a_plus, k1a_plus) instead of (a, k1a).
% nu_k1a at a_plus:
term1_prime_a = k1 * a_plus;
term2_prime_a = - exp(1i * pi/3) * xl * (k1*a_plus)^(1/3);
term3_prime_a = (1/60) * exp(1i * (2*pi/3)) * (xl^2) * ((1/2)*k1*a_plus)^(-1/3);
term4_prime_a = - ((xl^3 + 10) / 1400) * ((1/2)*k1*a_plus)^(-1);
term5_prime_a = - exp(1i * pi/3) * (281*xl^4 + 10440*xl) / 4536000 * ((1/2)*k1*a_plus)^(-5/3);
term6_prime_a = - exp(1i * 2*pi/3) * (73769*xl^5 + 6624900*xl^2) / 10478160000 * ((1/2)*k1*a_plus)^(-7/3);
term7_prime_a = ((93617*xl^6 + 16495400*xl^3 - 1744600) / 100900800000) * ((1/2)*k1*a_plus)^(-3);
term8_prime_a = (k1*a_plus)^(-11/3);
nu_k1a_plus = term1_prime_a + term2_prime_a + term3_prime_a + term4_prime_a + term5_prime_a + term6_prime_a + term7_prime_a + term8_prime_a;
% zeta_k1a at a_plus:
bracketTerm_plus = ...
    log( ( nu_k1a_plus + sqrt(nu_k1a_plus^2 - k1a_plus^2 ) ) / k1a_plus ) ...
    - sqrt( (nu_k1a_plus^2 - k1a_plus^2) / (nu_k1a_plus^2) );
zeta_k1a_plus = ( (3/2)^(2/3) ) * ( bracketTerm_plus )^(2/3);
% B0_zeta_k1a at a_plus:
B0_zeta_k1a_plus = -5/(48 * zeta_k1a_plus^2) ...
    + (1/(24*sqrt(zeta_k1a_plus))) * ...
    ( (5*nu_k1a_plus^3)/((nu_k1a_plus^2 - (k1*a_plus)^2)^(3/2)) ...
    - (3*nu_k1a_plus)/((nu_k1a_plus^2 - (k1*a_plus)^2)^(1/2)) );
% Hankel function at a_plus:
zArg_plus = exp(1i*(2*pi/3)) * (nu_k1a_plus^(2/3)) * zeta_k1a_plus;
AiVal_plus      = airy(0, zArg_plus);
AiPrimeVal_plus = airy(1, zArg_plus);
prefactor_plus = 2*sqrt(2) * ...
    ( (nu_k1a_plus^2 * zeta_k1a_plus) / (nu_k1a_plus^2 - (k1a_plus)^2 ) )^(1/4);
term1_br_plus = exp(-1i*pi/3)*AiVal_plus     * (nu_k1a_plus^(-1/3));
term2_br_plus = exp( 1i*pi/3)*AiPrimeVal_plus* (nu_k1a_plus^(-5/3)) * B0_zeta_k1a_plus;
bracket_plus  = term1_br_plus + term2_br_plus;
Term3_factor = 1 + 1/(nu_k1a_plus^2);   % if your O(1/nu^2) is based on the new nu_k1a_plus
H_plus = prefactor_plus * bracket_plus * Term3_factor;
% ==========================
% Evaluate at a_minus = a - da
% ==========================
a_minus   = a_orig - da;
k1a_minus = k1 * a_minus;
% --- Replicate the same lines, but with (a_minus, k1a_minus) ---
term1_prime_a_minus = k1 * a_minus;
term2_prime_a_minus = - exp(1i * pi/3) * xl * (k1*a_minus)^(1/3);
term3_prime_a_minus = (1/60) * exp(1i * (2*pi/3)) * (xl^2) * ((1/2)*k1*a_minus)^(-1/3);
term4_prime_a_minus = - ((xl^3 + 10) / 1400) * ((1/2)*k1*a_minus)^(-1);
term5_prime_a_minus = - exp(1i * pi/3) * (281*xl^4 + 10440*xl) / 4536000 * ((1/2)*k1*a_minus)^(-5/3);
term6_prime_a_minus = - exp(1i * 2*pi/3) * (73769*xl^5 + 6624900*xl^2) / 10478160000 * ((1/2)*k1*a_minus)^(-7/3);
term7_prime_a_minus = ((93617*xl^6 + 16495400*xl^3 - 1744600) / 100900800000) * ((1/2)*k1*a_minus)^(-3);
term8_prime_a_minus = (k1*a_minus)^(-11/3);
nu_k1a_minus = term1_prime_a_minus + term2_prime_a_minus + term3_prime_a_minus + term4_prime_a_minus + term5_prime_a_minus + term6_prime_a_minus + term7_prime_a_minus + term8_prime_a_minus;
bracketTerm_minus = ...
    log( ( nu_k1a_minus + sqrt(nu_k1a_minus^2 - k1a_minus^2 ) ) / k1a_minus ) ...
    - sqrt( (nu_k1a_minus^2 - k1a_minus^2) / (nu_k1a_minus^2) );
zeta_k1a_minus = ( (3/2)^(2/3) ) * ( bracketTerm_minus )^(2/3);
B0_zeta_k1a_minus = -5/(48 * zeta_k1a_minus^2) ...
    + (1/(24*sqrt(zeta_k1a_minus))) * ...
    ( (5*nu_k1a_minus^3)/((nu_k1a_minus^2 - (k1*a_minus)^2)^(3/2)) ...
    - (3*nu_k1a_minus)/((nu_k1a_minus^2 - (k1*a_minus)^2)^(1/2)) );
zArg_minus = exp(1i*(2*pi/3)) * (nu_k1a_minus^(2/3)) * zeta_k1a_minus;
AiVal_minus      = airy(0, zArg_minus);
AiPrimeVal_minus = airy(1, zArg_minus);
prefactor_minus = 2*sqrt(2) * ...
    ( (nu_k1a_minus^2 * zeta_k1a_minus) / (nu_k1a_minus^2 - (k1a_minus)^2 ) )^(1/4);
term1_br_minus = exp(-1i*pi/3)*AiVal_minus      * (nu_k1a_minus^(-1/3));
term2_br_minus = exp( 1i*pi/3)*AiPrimeVal_minus * (nu_k1a_minus^(-5/3)) * B0_zeta_k1a_minus;
bracket_minus  = term1_br_minus + term2_br_minus;
Term3_factor_minus = 1 + 1/(nu_k1a_minus^2);
H_minus = prefactor_minus * bracket_minus * Term3_factor_minus;
% ==========================
% Central difference
% ==========================
dH_da = (H_plus - H_minus) / (2*da);
% ==========================
% Display derivative
% ==========================
disp('Derivative of H_nu^(1)(k1*a) w.r.t. a is:');
disp(dH_da);

[Torsten]

The safest way is to write a function with input argument "a":

%% Constants
c = 3e8;                   % Speed of light
freq = 1e9;                % Frequency (1 GHz)
lambda0 = c / (freq);        % Free-space wavelength
k0 = 2 * pi ./ (lambda0);     % Free-space wavenumber
b = 20 ./ (k0);               % Outer radius of the dielectric coating (~0.955 meters)
d = 0.30*lambda0;% thickness of the dielectric material in cylider
a = b - d; % Compute inner radius of the cylinder
da = 1e-6;
Hprime = (fun(a+da)-fun(a-da))/(2*a)

Hprime = 4.8465e-08 - 1.0525e-07i

function H = fun(a)
%% Constants
mu0 = 4 * pi * 1e-7;       % Permeability of free space
eps0 = 8.854e-12;          % Permittivity of free space
c = 3e8;                   % Speed of light
freq = [1e9];                % Frequency (1 GHz)
lambda0 = c / (freq);        % Free-space wavelength
k0 = 2 * pi ./ (lambda0);     % Free-space wavenumber
b = 20 ./ (k0);               % Outer radius of the dielectric coating (~0.955 meters)
eps_r = 4;  
%mu_r= mu/mu0% Relative permittivity
eps1 = eps_r * eps0;       % Permittivity of dielectric
k1 = k0 * sqrt(eps_r);  
xl = -(2.33810741);
omega= 2*pi*freq;
d = 0.30*lambda0;% thickness of the dielectric material in cylider
%% Compute k1*a and k1*b
k1a = k1 * (a);
k1b = k1 * (b);
% 
%disp (k1a)
% %k1= omega* sqrt (eps1); %% if Naishadham's paper 
%% 
%% nu_equation for k1a case
term1 = k1 * a;
term2 = - exp(1i * pi/3) * xl * (k1*a)^(1/3);
% Term 3 
term3 = (1/60) ...
    * exp(1i * (2*pi/3)) ...
    * (xl^2) ...
    * ((1/2)*k1*a)^(-1/3);
% Term 4:
term4 = - ((xl^3 + 10) / 1400) ...
    * ((1/2)*k1*a)^(-1);
% Term 5:
term5 = - exp(1i * pi/3) ...
    * (281*xl^4 + 10440*xl) / 4536000 ...
    * ((1/2)*k1*a)^(-5/3);
% Term 6:
term6 = - exp(1i * 2*pi/3) ...
    * (73769*xl^5 + 6624900*xl^2) / 10478160000 ...
    * ((1/2)*k1*a)^(-7/3);
% Term 7:
term7 = ( (93617*xl^6 + 16495400*xl^3 - 1744600) / 100900800000 ) ...
    * ((1/2)*k1*a)^(-3);
% Term 8 (Big-O term):
term8 = (k1*a)^(-11/3);
nu_k1a = term1 + term2 + term3+term4+term5+term6+term7+term8;
%disp(nu_k1a);
%% write expression for zeta_k1a using equation 16
bracketTerm = ...
    log( ( nu_k1a + sqrt(nu_k1a^2 - k1a^2 ) ) / k1a ) ...
    - sqrt( (nu_k1a^2 - k1a^2) / (nu_k1a^2) );
zeta_k1a = ( (3/2)^(2/3) ) * ( bracketTerm )^(2/3);
%disp('zeta_k1a =');
%disp(zeta_k1a);
%% write beta0 zeta_k1a case
B0_zeta_k1a = -5/(48 * zeta_k1a^2) ...
    + (1 / (24 * sqrt(zeta_k1a))) * ...
    ( (5 * nu_k1a^3) / ((nu_k1a^2 - (k1*a)^2)^(3/2)) ...
    - (3 * nu_k1a)   / ((nu_k1a^2 - (k1*a)^2)^(1/2)) );
%disp('B0_zeta_k1a');
%disp(B0_zeta_k1a);
%% Hankel function first kind expansion k1a case 
% 2) Airy argument
% Now we have nu_k1a in place of nu, and zeta_k1a in place of zeta.
zArg = exp(1i * (2*pi/3)) * (nu_k1a^(2/3)) * zeta_k1a;
% 3) Compute Airy function and derivative
% airy(0,z) = Ai(z)
% airy(1,z) = Ai'(z)
AiVal      = airy(0, zArg);  % Ai(zArg)
AiPrimeVal = airy(1, zArg);  % Ai'(zArg)
% prefactor = 2*sqrt(2) * [ (nu_k1a^2 * zeta_k1a) / (nu_k1a^2 - (k1a)^2 ) ]^(1/4)
prefactor = 2*sqrt(2) * ...
    ( (nu_k1a^2 * zeta_k1a) / (nu_k1a^2 - (k1a)^2 ) )^(1/4);
% bracket = e^(-i*pi/3)*AiVal*nu_k1a^(-1/3)
%         + e^( i*pi/3)*AiPrimeVal*nu_k1a^(-5/3)*B0_zeta_k1a
term1 = exp(-1i * pi/3) ...
    * AiVal ...
    * (nu_k1a^(-1/3));
term2 = exp(1i * pi/3) ...
    * AiPrimeVal ...
    * (nu_k1a^(-5/3)) ...
    * B0_zeta_k1a;
bracket = term1 + term2;
Term3 = 1 + 1/(nu_k1a^2);
% ========================================
H_nu_k1a_approx = prefactor * bracket * Term3;
% ======================
% 8) Display the result
% ======================
%disp('Asymptotic Hankel H_nu^(1)(k1*a) = ');
%disp(H_nu_k1a_approx);
%% write the derivative of Hankel function first kind k1a case
%%%% [derivative of equation 21] 
%2) Define the Hankel function (first kind) as a function handle
%Hfun(a)
%2) Define the big-O replacement factor:
Term3 = 1 + 1/(nu_k1a^2);
% 3) Build a function handle for H_nu^(1)(k1*a)
%    We implement the prefactor, bracket, and Term3 as described above.
H = ( 2*sqrt(2) ...
    * ( (nu_k1a^2 * zeta_k1a) ./ (nu_k1a^2 - (k1.*a).^2 ) ).^(1/4) ) ...
    .* ( ...
    exp(-1i*pi/3)*AiVal* nu_k1a^(-1/3) ...
    + exp( 1i*pi/3)*AiPrimeVal*nu_k1a^(-5/3) * B0_zeta_k1a ...
    ).* Term3;
end

[SSWH2]

@Torsten thank you, I think this will work!

[SSWH2]

@Torsten can you please tell me how can Iwrite complex conjugate of this hankel function order 1 of my expression?[21 equation in screenshot] conj command will give direct value, but if I want to take derivative /prime of the conj command then it is showing wrror, so how can we write the complex conjugate?

[Torsten]

Where do you need to explicitly differentiate an expression in the computation of H ?

[Paul]

@Torsten

Shouldn't the denominator of this expression

Hprime = (fun(a+da)-fun(a-da))/(2*a)

be 2*da ?

[SSWH2]

Yes, you are right, I wrote my derivative like this-

dH_da = (H_plus - H_minus) / (2*da);

Now my question is for hankel functio second kind of the given expresssion , I need to do complex conjugate of this whole given equation, how can I do that?

[Torsten]

Sure, I didn't see this. I was already surprised about the small values for the derivative. Thank you @Paul .

[SSWH2]

@Torsten yup! now the problem is Conj command will not work for this case I think, beacuse I need second order hankel from the complex conjugate then also need derivative of the second order, ..don't know how to do .

[Sam Chak]

@SSWH2,

I am not overly familiar with the Hankel function, but I am curious. With a relatively large value for the real part of the complex order nu, is it reasonable for the Hankel functions of the first kind to produce a relatively small value?

%% Constants
c       = 3e8;          % Speed of light
freq    = 1e9;          % Frequency (1 GHz)
lambda0 = c/freq;       % Free-space wavelength
k0      = 2*pi/lambda0; % Free-space wavenumber
b       = 20/k0;        % Outer radius of the dielectric coating (~0.955 meters)
d       = 0.30*lambda0; % thickness of the dielectric material in cylider
a       = b - d         % Compute inner radius of the cylinder

a = 0.8649

[H1, nu_k1a] = fun(a)

H1 = -0.1275 + 0.2372i

nu_k1a = 40.0812 + 6.7302i

[Torsten]

@SSWH2

I don't understand what you exactly want. If f: IR -> C is a differentiable function, then d/dx ( conj(f) ) = conj ( d/dx(f) ). Does that help ?

[Sam Chak]

@SSWH2, Is the following a correct method to verify whether your code for the Hankel functions of the first kind returns the same value as HankelH1 = besselh(nu, 1, Z), when nu is a real constant?

%% Constants
c       = 3e8;          % Speed of light
freq    = 1e9;          % Frequency (1 GHz)
lambda0 = c/freq;       % Free-space wavelength
k0      = 2*pi/lambda0; % Free-space wavenumber
b       = 20/k0;        % Outer radius of the dielectric coating (~0.955 meters)
d       = 0.30*lambda0; % thickness of the dielectric material in cylider
a       = b - d         % Compute inner radius of the cylinder

a = 0.8649

[H1, nu_k1a] = fun(a)

H1 = -0.1338 - 0.0979i

nu_k1a = 2

HankelH1 = besselh(nu_k1a, 1, a)

HankelH1 = 0.0878 - 2.0744i

[SSWH2]

@Torsten fisrt of d/dx(f) is not working or diff(f) is not working in my code that's why I wrote manually the whole code as you can see in my code. Because its Hankel function but with complex nu , Conj command is working , first I need the complex conjugate of this Hankel function order 1, then I need the derivative of that complex conjugate function , which is basically complex number, so how can you take the derivative of a complex number in MATLAB? can you kindly suggest me?

[SSWH2]

@Sam Chak my code has no real number for nu, my nu is complex so it will not return same value with the built in besselh function of MATLAB.

[Torsten]

(conj(fun(a+da))-conj(fun(a-da)))/(2*da)

which should equal

conj(fun(a+da)-fun(a-da))/(2*da)

does not work ?

[SSWH2]

@Torsten the command you mentioned, will give me conjugate of the derivative, but I need complex conjugate of the Hankel function order 1 first, then I need the derivative of that complex conjugate , so what Iam thinking is that this Hankel function equation if you see my screenshot attached I have to make the complex conjugate manually by putting opoosite signs before i then I have to take the derivative as you suggested first time , isn't it ?

[Torsten]

The derivative of the complex conjugate of a function H: R -> C is equal to the complex conjugate of the derivative of the function. Thus you can compute the derivative of H first and then take the conjugate of the result.

syms a real
syms x
f = 2*x^2*a;
conj(diff(f,a))

ans = 

diff(conj(f),a)

ans = 

[SSWH2]

@Torsten do you mean taking complex conjugate of the derivative of the first kind (conj(fun(a+da))-conj(fun(a-da)))/(2*da)- will give me the derivative of the second kind using your command?

[Torsten]

I don't know what you mean by "derivative of the first kind" and "derivative of the second kind".

Do you mean first and second derivative of H with respect to a ? An approximation for the second derivative is given by

(fun(a+da)-2*fun(a)+fun(a-da))/da^2

[SSWH2]

@Torsten no its not just doing differentiation 2 times,, Hankel function first and second kind has different characteristics. I am telling, first kind of my given Hankel function expression become second kind when we will take the complex conjugate [ that's the rule], then how can I take the derivative of that complex equation using any built in command in matlab? Is there any built in funtion in Matlab to to take derivative of this?

[Torsten]

Here are the relations for the derivative of the Hankel-functions of the first and of the second kind (expressed by the Hankel-functions themselves):

https://mathworld.wolfram.com/HankelFunctionoftheFirstKind.html

https://mathworld.wolfram.com/HankelFunctionoftheSecondKind.html

[SSWH2]

@Torsten Thanks for your prompt reply! But this will not work my equation nu /n (order) of the Hankel function is not real its complex, that's why I am looking for any built in function to do complex conjugate and derivative.

[Torsten]

Then I suggest you compute the Bessel functions Jn and Yn separately and form the Hankel functions of first and second kind afterwords.

Here is code to compute the Bessel functions Jn and Yn of complex order and argument:

https://uk.mathworks.com/matlabcentral/fileexchange/9515-bessel-function-of-complex-order-and-argument

[SSWH2]

My equations are already given to Hankel, conversion to bessel and deffirentitate will not work in this case , I tried with the given complex order code long time ago, this will not work I thnk, anyways thanks for the reply !

[Torsten]

Why ? Differentiate J and Y you receive from the code I gave you the link for by a finite difference quotient and compute the Hankel derivative of first and second kind as dJ + 1i*dY resp. dJ - 1i*dY.

[SSWH2]

@Torstenmanuel derivative I already did, so I thought is there any buit in function in MATLAB which can does it randomly. I did already manually this way taking small difference at a and b. Thanks

% Central difference

% 1) Save your original b, k1b

b_orig = b;

k1b_orig = k1b; % presumably k1*b

% 2) Choose a small finite-difference step

db = 1e-6;

% ==========================

b_minus = b_orig - db;

dH_db = (H_plus_k1b - H_minus_k1b) / (2 * db);