A quatity is being solved by a self consistent integration

How to find Z from the below equation

here \mathcal{P} means the principal value integration. I was tried in the following way, but couldn't figure out how to solve this,

A = 2000; a = 500; tolerance = 10^-4; Z = 0;
for i = 1 : 10
    result = integral(@(x) (x.^2.*((A^2+Z^2)./(A^2+((x.^2+a^2)))) .* (sqrt(x.^2+a^2).*(x.^2+a^2-Z^2)).^(-1)), 0,A, 'PrincipalValue', true);
    new_Z = sqrt(result);
    if abs(new_Z - Z) < tolerance
        Z = new_Z;
        break;
    end
    Z = new_Z;
end
disp(new_Z);

Thank you in advance!

[Torsten]

Why is the Principal Value necessary to be taken ? In case a^2 - Z^2 <= 0 ?

[pritha]

Hi Torsten,

Z is completely unknown and there was no specified situation for that.

Answer 1

[Torsten]

format long
syms x 
A = 2000;
a = 500;
b = 1000;
Z = 0;
for i=1:20
  f = x^4*((A^2+Z^2)/(A^2+4*(x^2+a^2)))^4 / (sqrt(x^2+a^2)*(x^2+a^2-Z^2));
  I = double(int(f,x,0,A,'PrincipalValue',true));
  Zpi = sqrt(b^2-I)
  Z = real(Zpi)
end

Zpi = 
     9.822515593894991e+02

Z = 
     9.822515593894991e+02

Zpi = 
      9.926239912427430e+02 -1.331193739501012e-147i

Z = 
     9.926239912427430e+02

Zpi = 
      9.942989166123056e+02 -4.915924275823428e-153i

Z = 
     9.942989166123056e+02

Zpi = 
      9.945686435588058e+02 +9.829182155008009e-152i

Z = 
     9.945686435588058e+02

Zpi = 
      9.946120599497613e+02 -1.207757180393064e-148i

Z = 
     9.946120599497613e+02

Zpi = 
      9.946190479156575e+02 -2.457171010268977e-153i

Z = 
     9.946190479156575e+02

Zpi = 
      9.946201726310163e+02 +1.228584115851398e-152i

Z = 
     9.946201726310163e+02

Zpi = 
      9.946203536539656e+02 +9.828671137972504e-153i

Z = 
     9.946203536539656e+02

Zpi = 
      9.946203827896022e+02 -2.061221673054303e-146i

Z = 
     9.946203827896022e+02

Zpi = 
      9.946203874789816e+02 +1.179440496446327e-151i

Z = 
     9.946203874789816e+02

Zpi = 
      9.946203882337369e+02 -1.106609953389711e-137i

Z = 
     9.946203882337369e+02

Zpi = 
      9.946203883552148e+02 -3.542724730738004e-147i

Z = 
     9.946203883552148e+02

Zpi = 
      9.946203883747665e+02 +1.006455889394421e-149i

Z = 
     9.946203883747665e+02

Zpi = 
     9.946203883779135e+02

Z = 
     9.946203883779135e+02

Zpi = 
      9.946203883784200e+02 -1.610329423025159e-147i

Z = 
     9.946203883784200e+02

Zpi = 
      9.946203883785015e+02 +1.228583849353811e-152i

Z = 
     9.946203883785015e+02

Zpi = 
      9.946203883785146e+02 -1.030610830736004e-145i

Z = 
     9.946203883785146e+02

Zpi = 
      9.946203883785167e+02 -2.061221661472003e-146i

Z = 
     9.946203883785167e+02

Zpi = 
      9.946203883785171e+02 +4.221381962694661e-143i

Z = 
     9.946203883785171e+02

Zpi = 
      9.946203883785171e+02 +1.376013911276247e-151i

Z = 
     9.946203883785171e+02

f = x^4*((A^2+Z^2)/(A^2+4*(x^2+a^2)))^4 / (sqrt(x^2+a^2)*(x^2+a^2-Z^2));
double(Z^2 - b^2 + real(int(f,x,0,A,'PrincipalValue',true)))

ans = 
    -6.035923459074367e-11

[pritha]

Hi Torsten,

Thank you. This works very well. However, when I run the code for 1500 values of 'a', it takes too much time, almost like 1hr. Could you please suggest me some wayout or any otherr process with which such kind of problem can be solved?

[Torsten]

If the values for A don't change much, you should use the result for Z of the call for A(i) as initial guess for the call with A(i+1).

Further, you could try to solve your equation directly without fixed-point iteration using the "vpasolve" function:

syms Z x 
A = 2000;
a = 500;
b = 1000;
f = x^4*((A^2+Z^2)/(A^2+4*(x^2+a^2)))^4 / (sqrt(x^2+a^2)*(x^2+a^2-Z^2));
eqn = Z^2 - b^2 + real(int(f,x,0,A,'PrincipalValue',true)) == 0;
vpasolve(eqn,Z)