how to solve transcedental equation in matlab

I am practicing solving the next transcendental equation in matlab

(a/c)*sqrt((b*m1)^2-(p*c)^2)-atan( sqrt(( (p*c)^2-(b*m2)^2  ) /( (b*m1)^2-(p*c)^2  ) ) )-atan( sqrt(( (p*c)^2-(b*m3)^2  ) /( (b*m1)^2-(p*c)^2  ) ) ) == r*pi

here

a=1x10^-6;
c= 3x10^8;
m1=2.2;
m2=1.5;
m3=1;

I was trying to plot "p" vs "b" where "b" runs from 0 to 3x10^15 and r is a parameter that takes values of 0, 1 and 2. I already tried all day but I cannot find solution, I tried with fzero(fun,xo) without success, can you give any suggestion?

[Walter Roberson]

My tests suggest that for r = 0, that the smallest roots are at b = +/- 3000000000000000*sqrt(259)*arctan(5*sqrt(259)*sqrt(5)*(1/259))*(1/259) and p = +/- (b/200000000), for a total of two roots for that b.

It appears that there might perhaps be one primary root for each b larger than that; I have not chased to see if the negative is also a root. I do not have a characterization of the roots for larger b. At most p about 2.182E7 for b = 3E15. Possibly approximately p = sqrt(b/(2*pi)) being slightly too high.

[Walter Roberson]

You ask about plotting p vs b where b runs from 0 to 3E15, but the axis that runs to about 3E15 on the plot is the left hand axis. That suggest that the plot is b vs p instead of p vs p.

It is possible to demonstrate that if a root exists at all, that it is in the range p = b/200000000 to 11*b/1500000000 so you could use those as the bounds on your fzero.

[alburary daniel]

without success today, Can you help me?

[Walter Roberson]

Please do not close questions that have an answer. You were informed about that before https://www.mathworks.com/matlabcentral/answers/405170-how-to-solve-transcendental-equations-in-matlab#comment_578186

I spent several hours working on your question, and it is very frustrating that my answer just disappeared.

Answer 1

[Walter Roberson]

I got code to work... but the trend is exactly opposite of what you are looking for.

Q = @(v) sym(v, 'r');
arctan = @atan;

a = Q(1*10^-6);
c = Q(3*10^8);
m1 = Q(2.2);
m2 = Q(1.5);
m3 = Q(1);
syms b p r
eqn = (a/c)*sqrt((b*m1)^2-(p*c)^2)-atan( sqrt(( (p*c)^2-(b*m2)^2  ) /( (b*m1)^2-(p*c)^2  ) ))-atan( sqrt(( (p*c)^2-(b*m3)^2  ) /( (b*m1)^2-(p*c)^2  ) ) ) - r*pi;

B0min = Q(3000000000000000)*sqrt(Q(259))*arctan(5*sqrt(Q(259))*sqrt(Q(5))*(1/Q(259)))*(1/Q(259));
B0max = Q(3*10^15);
B0 = linspace(B0min, B0max, 100);

nB0 = length(B0);
R = 0 : 2;
nR = length(R);

ps = zeros(nB0, nR, 'sym');

for ridx = 1 : nR
    this_r = R(ridx);
  eqnr = subs(eqn, r, this_r);
    for K = 1 : nB0
        this_b = B0(K);
        this_eqn = subs(eqnr, b, this_b);

          sols = vpasolve(this_eqn, p, [this_b/200000000, 11*this_b/1500000000]);

          if isempty(sols)
              fprintf('No symbolic solution for eqn, r = %d, b = %g\n', this_r, this_b);
              sols = nan;
          else
              % fprintf('symbolic okay for eqn, r = %d, b = %g\n', this_r, this_b);
          end
          ps(K, ridx) = sols;
      end
  end

plot(ps, B0)
xlabel('p')
ylabel('b')
legend( sprintfc('r = %d', R) )

[alburary daniel]

but what happens if the necesary condition to find the roots is that m2*b<p<m1*b for any value of b?

[Walter Roberson]

"but what happens if the necesary condition to find the roots is that m2*b<p<m1*b for any value of b?"

Then there are no roots. I proved mathematically that the condition for real valued roots is

200000000*b <= p <= 1500000000*b

and the negatives of those values.

Values in the range 1.5*b to 2.2*b are for certain not possible solutions.

[alburary daniel]

So, this condition that you mention before it could be changed, for example, if now m1=1.2 and m2=m3=1, or remains equal?

[Walter Roberson]

The lower limit is m2/300000000 . The upper limit is m1/300000000

More precisely, the lower limit is m2/c and the upper limit is m1/c .

If you are expecting values in the range of m1*b to m2*b then your c would need to be approximately 1.

[alburary daniel]

so the lower limit is B0min = Q(3000000000000000)*sqrt(Q(259))*arctan(5*sqrt(Q(259))*sqrt(Q(5))*(1/Q(259)))*(1/Q(259)); for the first one?

[Walter Roberson]

B0min = arctan(sqrt((m2-m3)*(m2+m3)/((m1-m2)*(m1+m2))))*c/(a*sqrt((m1-m2)*(m1+m2)))

[Walter Roberson]

Here is an updated version of the code:

Q = @(v) sym(v, 'r');
arctan = @atan;

a = Q(1*10^-6);
c = Q(3*10^8);
m1 = Q(2.2);
m2 = Q(1.5);
m3 = Q(1);

syms b p r
eqn = (a/c)*sqrt((b*m1)^2-(p*c)^2)-atan( sqrt(( (p*c)^2-(b*m2)^2  ) /( (b*m1)^2-(p*c)^2  ) ))-atan( sqrt(( (p*c)^2-(b*m3)^2  ) /( (b*m1)^2-(p*c)^2  ) ) ) - r*pi;

%the following B0min is for 2.2, 1.5, 1
%B0min = Q(3000000000000000)*sqrt(Q(259))*arctan(5*sqrt(Q(259))*sqrt(Q(5))*(1/Q(259)))*(1/Q(259));
B0min = arctan(sqrt((m2-m3)*(m2+m3)/((m1-m2)*(m1+m2))))*c/(a*sqrt((m1-m2)*(m1+m2)));

B0max = Q(3*10^15);
B0 = linspace(B0min, B0max, 100);

nB0 = length(B0);
R = 0 : 2;
nR = length(R);

ps = zeros(nB0, nR, 'sym');

for ridx = 1 : nR
    this_r = R(ridx);
    eqnr = subs(eqn, r, this_r);
    badrun = false;
    for K = 1 : nB0
        this_b = B0(K);
        this_eqn = subs(eqnr, b, this_b);
        sols = vpasolve(this_eqn, p, [this_b*m2/c, this_b*m1/c]);
        % "the necesary condition to find the roots is that m2*b<p<m1*b for any value of b"
        %sols = vpasolve(this_eqn, p, [this_b*m2, this_b*m1]);

          if isempty(sols)
              if ~badrun
                  fprintf('No symbolic solution for eqn starting from r = %d, b = %g\n', this_r, this_b);
                  badrun = true;
              end
              sols = nan;
          else
              if badrun
                  fprintf('symbolic okay again starting from r = %d, b = %g\n', this_r, this_b);
              elseif K == 1
                  fprintf('symbolic okay starting from r = %d, b = %g\n', this_r, this_b);

              end
              badrun = false;
          end
          ps(K, ridx) = sols;
      end
      if badrun
          fprintf('No symbolic solution for eqn through to r = %d, b = %g\n', this_r, this_b);
          badrun = false;
      else
          fprintf('symbolic okay through to r = %d, b = %g\n', this_r, this_b);
      end
  end

plot(ps, B0)
xlabel('p')
ylabel('b')
legend( sprintfc('r = %d', R) )

[alburary daniel]

I´m still a doubt about what happens in the case when m2=m3 and r=0, the minimum is zero, rigth?

[Walter Roberson]

The case where m2=m3 and r=0 has a theoretical solution of b=0, p=m2/c*b = m2/c*0 = 0 . This requires taking the limiting case, as substituting b=0, p=0 directly will get you the indeterminate 0/0

[alburary daniel]

but I don't understand if above condition is indeterminate how can made this type of graph?

[Walter Roberson]

It is not clear how the variable names there relate to your m1, m2, m3 ?

[alburary daniel]

well, here n1=m1 and n2=m2=m3

[Walter Roberson]

I do not know. I notice that they are using different c values in the diagram, but you set up your equation with only one c value.