Using Matlab to find numerically all roots and plotting them
5 views (last 30 days)
Show older comments
I want to solve a trigonometric equation numerically for a via Matlab
sin(k+a)=xcos(k)
I want to have solution sol[x,k], so that for each x(specially with this one) and k I get (either zero, one, two or more) solutions a. Then to plot them, or visualize via a vector or table. Table of values to explicitly see the solutions, i.e., which ranges of x and k had one, two, three or no solutions.
The code that I wrote upto now is
clear all
k=[0:0.1:pi]
x=[0:0.1:pi]
f=sin(k+a)-x*cos(k);
for i=1:3
vpasolve(f,a,'random',true)
end
0 Comments
Answers (1)
John D'Errico
on 9 Dec 2018
Edited: John D'Errico
on 9 Dec 2018
There are infinitely many solutions. So wanting to find all of them might take some time. But solve can easily find some principal solutions.
syms a k x
a = solve(sin(k+a) == x*cos(k),a)
a =
asin(x*cos(k)) - k
pi - k - asin(x*cos(k))
Pencil and paper would suffice as easily, since you can just take the inverse sine of both sides. I'm not at all sure why you want to go through the effort that you did.
Just add any integer multiple of 2*pi to find all solutions in general, since the sine function is periodic with period 2*pi.
2 Comments
John D'Errico
on 9 Dec 2018
Edited: John D'Errico
on 9 Dec 2018
As you can see from the analytical solution, there will in general be exactly two principal solutions. Rarely, we will see exactly one solution. (Or no real solutions, for some values of x.) Again, given any solution, you can add n*2*pi to it to gain another solution, but any integer value of n.
Now, can you find exactly 3 solutions in a closed interval of length 2*pi? Well, yes, in case one of the solutions lies eactly on one end of the interval.
But, yes, you can have a scenario where there are no real solutinos at all.
For example, choose k==0, and x==2. Then x*cos(k) will be 2. What are the real solutions of
sin(x+a) == 2
There are no real solutions. Period.
Again, the simple solution is to just use what solve gives. And there is absolutely no reason to use a loop, vpasolve, etc.
k = 0:0.1:pi;
x = (0:0.4:pi)';
a1 = asin(x.*cos(k)) - k;
a2 = pi - k - asin(x.*cos(k));
So a1 and a2 are the two distinct solutions. We can exclude the complex solutions simply enough, by just setting them to NaN.
wascomplex = imag(a1) ~= 0;
a1(wascomplex) = NaN;
a2(wascomplex) = NaN;
Here are the TWO solutions, as a function of x and k:
surf(k,x,a1)
hold on
surf(k,x,a2)
xlabel 'k'
ylabel 'x'
zlabel 'a1 & a2'

The one scenario where you will find 3 solutions in a closed interval of length 2*pi? See what happens when x=k=0.
Then we have solutions for a at exactly [0, pi, 2*pi]. But that is only because we can freely add 2*pi to the solution at 0 to get another solution.
Oh, under what scanario will you have only one solution? For that to happen, we need to have a1 == a2, in the above solution. then...
asin(x.*cos(k)) - k = pi - k - asin(x.*cos(k))
a little quick algebra will reduce that to
2*asin(x.*cos(k)) == pi
Divide by 2, then take the sin of both sides, recognizing that sin(pi/2) = 1. So those two solutions reduce to one only when we have
x*cos(k) == 1
If abs(x) is even slightly greater in magnitude than abs(1/cos(k)), even by eps, then the solution goes complex, so no real solutions.
Using a little mathematics is always a better choice than just throwing a problem at vpasolve, abdicating all thought.
See Also
Categories
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!