Numerically solving x + Asin(Wx+C) = 0, for x. (A, W & C are known constants)
2 views (last 30 days)
Show older comments
I need to solve for x in "x + Asin(Wx+C) = 0", where A, W and C are all known constants. I looked and this seems to have no analytical solution. Hence, I tried Newton-Rapson Method. I managed to achieve convergence by rearranging the equation in the form A/x + 1/sin(Wx+c). However, when I substitute back for the answer for some instances there is a significant error. In order to recreate this seemingly random error I change, C and perform the evaluation. I suspect this is due to the numerical precision issues. I've taken some steps to mitigate the numerical precision issues. How can I reduce the error?
%Numerically Solving x + A*sin(W*x + C);
%Solving 0 = A/x + 1/sin(W*x + C) for fast convergence
close all; clear; clc;
%Specification
A = 1.6E4;
W = 1.57E-9;
C0 = 1.5E-5;
Two_Pi = 2*pi;
%Specifying the Numerical Uncertainty
NUN = 1E-18; %In Calculation
%Simulation Setup
Ns = 1E6;
N0 = 9.1433E9;
%Initializing
%The Error Vector
Err_Est_Vec = zeros(Ns,1);
%Iteration Counter
Itt_Count_Vec = zeros(Ns,1);
for n = 1 : Ns
C = mod(C0*(N0+n),Two_Pi);
%Initial Iteration Solution Estimate
x_C = -A*sin(C);
%Iteration Counter
Itt_Count = 0;
while (abs(x_C + A*sin(mod(W*x_C + C,Two_Pi))) > NUN) && (Itt_Count < 1E3)
Itt_Count = Itt_Count + 1;
x_N = x_C + x_C*sin(mod(W*x_C+C,Two_Pi))*(x_C + A*sin(mod(W*x_C+C,Two_Pi)))...
./(A*(sin(mod(W*x_C+C,Two_Pi))).^2 + W*cos(mod(W*x_C+C,Two_Pi))*(x_C.^2));
%Updating the Next Estimate as the Current
x_C = x_N;
end %WHILE
%Updating Iteration Counter Vector & Error Estimate
Itt_Count_Vec(n) = Itt_Count;
Err_Est_Vec(n) = x_C + A*sin(mod(W*x_C + C,Two_Pi));
end %FOR n
figure;
subplot(211); stem(Err_Est_Vec); title('Error');
subplot(212); stem(Itt_Count_Vec); title('Iteration Count');
0 Comments
Answers (3)
Star Strider
on 14 Nov 2018
2 Comments
Star Strider
on 14 Nov 2018
My pleasure.
I do not know what options you have available. A Symbolic Math Toolbox solution would be:
syms x
A = 1.6E4;
W = 1.57E-9;
C0 = 1.5E-5;
Eq = x + A*sin(W*x+C0) == 0;
X = vpasolve(Eq,x)
producing:
X =
-0.23999397134244056217948413161828
John D'Errico
on 14 Nov 2018
Edited: John D'Errico
on 14 Nov 2018
When you formulate the problem by inversion as you did, you might expect to see numerical problems some of the time, depending on the value of x. That is not even a remote surprise to me.
x + Asin(Wx+C) = 0
Look at the constants involved.
A = 1.6E4;
W = 1.57E-9;
C0 = 1.5E-5;
Large numbers and small numbers are an assurance of nastiness in a numerical problem when you are working in floating point arithmetic. Then inverting things is just asking for problems. How would I solve this? The obvious answer is to use the symbolic toolbox.
vpasolve(x + A*sin(W*x+C0),x)
ans =
-0.23999397134244056217948413161828
FZERO has no problems in finding that solution.
fzero(@(x) x + A*sin(W*x+C0),0)
ans =
-0.239993971342441
Are you looking for the smallest positive solution? There don't seem to be any that I see. Lets see, if I differentiate that expression, the function has a monotonic increasing first derivative for all positive x.
vpa(diff(x + A*sin(W*x+C0)),5)
ans =
0.00002512*cos(1.57e-9*x + 0.000015) + 1.0
So for all positive x, there is provably no zero.
Is there a simple solution? There very much is. The clue comes from your constants. So try this:
T = taylor(x + A*sin(W*x+C0),'expansionpoint',C0);
vpa(T,3)
ans =
1.0*x - 2.96e-19*(x - 1.5e-5)^2 - 1.03e-23*(x - 1.5e-5)^3 + 6.08e-38*(x - 1.5e-5)^4 + 1.27e-42*(x - 1.5e-5)^5 + 0.24
Think about what I did. I expanded the problem as a Taylor series around C0. C0 dominates the argument of sin there, therefore we should expect that series to be extremely rapidly convergent. In fact, the problem now becomes trivial.
You can recognize that ALL terms past the first order term are completely insignificant. TOSS THEM!!!!!!! What does the problem reduce to now? Yep. In 16 significant digits, we get:
1.000025119999997*x + 0.239999999991 == 0
Can you solve that problem. Do you really need fzero? ;-)
x = -0.239999999991/1.000025119999997
x =
-0.239993971342441
That looks just like what we got from vpasolve.
This does not mean you should be just using vpasolve. But an intelligent use of a series expansion seems very much to the point, and solves the problem trivially.
0 Comments
Image Analyst
on 15 Nov 2018
Edited: Image Analyst
on 15 Nov 2018
One very simple way to solve for x when y is zero is to just create a bunch of y and see which is closest to zero. No fancy, complicated functions just a simple formula followed by min(abs(y)). Basically 3 lines of code: one to setup x, one to create y from x and the constants, and one to find the zero crossing. Here is the full demo:
numberOfPoints = 1000000;
% Smaller numberOfPoints for less precision.
% Larger numberOfPoints for more precision.
A = 1.6E4;
W = 1.57E-9;
C = 1.5E-5;
% 2 * pi / period = W, so 2 * pi/W = period.
period = 2 * pi/W % = 4002028858.07617 !!!
% Plot function over 4 periods.
x = linspace(-2*period, 2*period, numberOfPoints);
y = x + A * sin(W * x + C);
subplot(2, 1, 1);
plot(x, y, 'b-', 'LineWidth', 2);
grid on;
fontSize = 18;
title('Zoomed out view of y = x + A * sin(W * x + C)', 'FontSize', fontSize);
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
% Hmmm.... looks like a straight line because A is so much smaller than x.
% Plot over a range closer to where the signal equals zero:
x = linspace(-0.26, -0.21, numberOfPoints);
y = x + A * sin(W * x + C);
subplot(2, 1, 2);
plot(x, y, 'b-', 'LineWidth', 2);
grid on;
title('Zoomed in view of y = x + A * sin(W * x + C)', 'FontSize', fontSize);
xlabel('x', 'FontSize', fontSize);
ylabel('y', 'FontSize', fontSize);
hold on;
ax = gca;
ax.XAxisLocation = 'origin';
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0, 0.04, 1, 0.96]);
% Find out where the equation is closest to zero
[minValue, indexAtMin] = min(abs(y))
message = sprintf('y is closest to zero at x = %f where y has a value of %e.', ...
x(indexAtMin), minValue);
fprintf('%s\n', message);
text(x(indexAtMin), -0.01, message, 'Color', 'r', 'FontSize', 16);
uiwait(helpdlg(message));
Here is the result:
![0000 Screenshot.png](https://www.mathworks.com/matlabcentral/answers/uploaded_files/195082/0000%20Screenshot.png)
The program says
y is closest to zero at x = -0.239994 where y has a value of 8.651757e-09.
which is the same result that John got.
0 Comments
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!