## how to solve transcedental equation in matlab

Asked by alburary daniel

### alburary daniel (view profile)

on 13 Jun 2018
Latest activity Edited by Walter Roberson

### Walter Roberson (view profile)

on 27 Jul 2018
Accepted Answer by Walter Roberson

### Walter Roberson (view profile)

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?

Show 1 older comment
Walter Roberson

### Walter Roberson (view profile)

on 13 Jun 2018
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

### alburary daniel (view profile)

on 15 Jun 2018
without success today, Can you help me?
Walter Roberson

### Walter Roberson (view profile)

on 27 Jul 2018
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.

### Tags ### Walter Roberson (view profile)

Answer by Walter Roberson

### Walter Roberson (view profile)

on 15 Jun 2018

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) )

Walter Roberson

### Walter Roberson (view profile)

on 22 Jun 2018
It is not clear how the variable names there relate to your m1, m2, m3 ?
alburary daniel

### alburary daniel (view profile)

on 23 Jun 2018
well, here n1=m1 and n2=m2=m3
Walter Roberson

### Walter Roberson (view profile)

on 24 Jun 2018
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.