Numerically solving x + Asin(Wx+C) = 0, for x. (A, W & C are known constants)

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

Answer 1

[Star Strider]

I am not certain what you are doing.

I would use fzero or fsolve.

[Thushara]

Thanks ST!

I didn't know about these functions. Thus, I did some research and got familier with 'fzero'. Unfortunately, 'fsolve' needs the Optimization Toolbox, which I don't have. With regards to 'fzero' my breif reserch show that it utilizes a number of method including the 'secant' method a varient of the 'Newton-Raphson' method I used in my example. It was quite interesting as both 'fzero' and my method failed in the same instances.

Regards,

T

[Star Strider]

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

Answer 2

[John D'Errico]

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.

Answer 3

[Image Analyst]

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:

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.