difference between simulating a polynomial to the nth power and using the polynomial expansion theorem to expand it. the index increases, the difference in values increases

%% Take the 2nd power expansion as an example
r_min=500000;
r_max=13242000;
gama=2.554325099941455e-16;
x= @(r0) sqrt(2).*0.1.*r0.^2; 
r0_values = linspace(r_min, r_max, 1000); 
a=arrayfun(@(r0) (1-2*exp(-x(r0).*gama)+exp(-2*x(r0).*gama))^2,r0_values) %

a = 1×1000

1.0e-08 * 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

b=arrayfun(@(r0) 1+4*exp(-2*x(r0).*gama)+exp(-4*x(r0).*gama)+ ...    %There are 6 items in the expanded type
    2*exp(-2*x(r0).*gama)-4*exp(-x(r0).*gama)-4*exp(-3*x(r0).*gama),r0_values)

b = 1×1000

1.0e-08 * 0 -0.0000 -0.0000 -0.0000 0 0.0000 -0.0000 0.0000 -0.0000 0 -0.0000 -0.0000 0.0000 0.0000 0.0000 -0.0000 -0.0000 0 -0.0000 -0.0000 0.0000 -0.0000 -0.0000 0.0000 0.0000 -0.0000 -0.0000 0 0 -0.0000

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

max(abs(a-b))

ans = 1.7406e-15

[Torsten]

I don't understand what you are asking for.

You mean you get different results depending on whether you use the expanded or the unexpanded form ? I wouldn't be surprised with such extreme values for gama, r_min and r_max.

[raven]

Yes, and the difference between them is significant as the number of times increases (not as close as to the index of 2).

The extrema of the parameters are derived from my research, do you mean to say that this is the reason for the divergence between the two?

My research is more complex, and these values need to be followed by integration and other operations, and the results show that if you the expanded’s results will converge as the exponent increases; And the unexpanded one is not

Is there any way to reduce this error?

[Torsten]

For the example given, the maximum difference between a and b is 1e-15 (see above).

Do you have an example where the difference is significant ?

And why do you want to use the expanded form ? Intuitively, its precision will always be worse than for the compact form.

[raven]

I have to use the expansion method in my research.

a and b need to be further multiplied by a same function, and then integrals with respect to r0,

In summary, the final results of both from the power 0 to the power of 30 are shown below, and it can be seen that the former converges with the increase of the number of powers, while the latter does not

[Torsten]

I still don't get why you need to expand the expression.

Answer 1

[Steven Lord]

Let's look at the numbers involved in the calculations in your polynomial evaluations. I've renamed r0_values to just r0 to make the expressions used to compute those numbers closer to the actual code used to create a and b. I also left out the constant 1 in termsInB and in the Term* variables in the results table, but added it in when I computed B in the results table.

%% Take the 2nd power expansion as an example
r_min=500000;
r_max=13242000;
gama=2.554325099941455e-16;
x= @(r0) sqrt(2).*0.1.*r0.^2; 
r0 = linspace(r_min, r_max, 1000).'; 
valuesToBeSquared = 1-2*exp(-x(r0).*gama)+exp(-2*x(r0).*gama);
termsInB = [4*exp(-2*x(r0).*gama), exp(-4*x(r0).*gama), ... 
            2*exp(-2*x(r0).*gama), -4*exp(-x(r0).*gama), -4*exp(-3*x(r0).*gama)];
tableB = array2table(termsInB, VariableNames = "Term" + (1:5));
results = [table(r0, valuesToBeSquared), tableB];
results = addvars(results, results.valuesToBeSquared.^2, ...
    'Before', 'valuesToBeSquared', 'NewVariableNames', 'A');
results = addvars(results, 1+sum(termsInB, 2), 'After', 'A', 'NewVariableNames', 'B');
format short
head(results)

        r0            A              B         valuesToBeSquared    Term1      Term2     Term3      Term4      Term5 
    __________    __________    ___________    _________________    ______    _______    ______    _______    _______

         5e+05    6.6515e-21              0       8.1557e-11        3.9999    0.99996         2         -4    -3.9999
    5.1275e+05    8.1364e-21    -4.4409e-16       9.0202e-11        3.9999    0.99996         2         -4    -3.9999
    5.2551e+05    9.9038e-21    -4.4409e-16       9.9518e-11        3.9999    0.99996         2         -4    -3.9999
    5.3826e+05    1.1998e-20    -4.4409e-16       1.0954e-10        3.9999    0.99996         2         -4    -3.9999
    5.5102e+05    1.4471e-20              0       1.2029e-10        3.9999    0.99996         2         -4    -3.9999
    5.6377e+05    1.7378e-20     4.4409e-16       1.3182e-10        3.9999    0.99995         2         -4    -3.9999
    5.7653e+05    2.0784e-20    -4.4409e-16       1.4417e-10        3.9999    0.99995         2         -4    -3.9999
    5.8928e+05     2.476e-20     4.4409e-16       1.5735e-10        3.9999    0.99995    1.9999    -3.9999    -3.9998

The values in the Term* variables smacks of catastrophic cancellation, or as MathWorks chief scientist Cleve Moler called it in this article from 1996 "severe subtractive cancellation." While the values in the Term* variables are not that large in an absolute sense, compared to the values in B they are relatively very large.

[raven]

I get it !

Thank you sooo much!

But, is there any way to reduce " catastrophic cancellation" ?

Although A and B are numerically close, since they need to be integrated regard to r0 later, when the exponent is large, the final result is very different

Answer 2

[Umar]

Hi @Raven,

After going through your comments and analysis of your code, there are two distinct approaches: simulating a polynomial raised to the nth power and using the polynomial expansion theorem. Each method has its own implications, especially as the degree of the polynomial increases. When you simulate a polynomial to the nth power, you directly compute the value of the polynomial for a given input. This approach involves evaluating the polynomial expression as it is defined, which can be computationally intensive for high degrees. The simulation yields precise results for specific values of the input variable. On the other hand, the polynomial expansion theorem allows us to express a polynomial in a more manageable form, often as a series expansion. This method can simplify calculations, especially for higher degrees, by breaking down the polynomial into a sum of terms. The expansion theorem can provide insights into the behavior of the polynomial as the degree increases. In the context of your question, as the index increases, the difference in values between the two methods becomes more pronounced. This is due to the nature of polynomial growth and the approximation errors that can arise from using series expansions. Let me illustrate this with the provided MATLAB code, which simulates the second power expansion and compares it with the polynomial expansion theorem.

% Define parameters
r_min = 500000;
r_max = 13242000;
gamma = 2.554325099941455e-16;

% Define the function x(r0)
x = @(r0) sqrt(2) .* 0.1 .* r0.^2;

% Generate r0 values
r0_values = linspace(r_min, r_max, 1000);

% Simulate the polynomial to the 2nd power
a = arrayfun(@(r0) (1 - 2 * exp(-x(r0) .* gamma) + exp(-2 * x(r0) .* gamma))^2, 
r0_values);

% Use polynomial expansion theorem
b = arrayfun(@(r0) 1 + 4 * exp(-2 * x(r0) .* gamma) + exp(-4 * x(r0) .* gamma)   
+ ... 2 * exp(-2 * x(r0) .* gamma) - 4 * exp(-x(r0) .* gamma) - 4 * exp(-3 *   
 x(r0) .* gamma), r0_values);

% Plotting the results
figure;
plot(r0_values, a, 'b-', 'DisplayName', 'Simulated Polynomial (2nd Power)');
hold on;
plot(r0_values, b, 'r--', 'DisplayName', 'Polynomial Expansion');
xlabel('r0 Values');
ylabel('Polynomial Value');
title('Comparison of Polynomial Simulation and Expansion');
legend show;
grid on;

% Display final results
disp('Final Results:');
disp('Simulated Polynomial Values (a):');
disp(a);
disp('Expanded Polynomial Values (b):');
disp(b);

Please see attached.

Again, the choice between simulating a polynomial to the nth power and using the polynomial expansion theorem depends on the specific requirements of your analysis. While direct simulation provides exact values, the expansion theorem can offer a more efficient approach for higher degrees, albeit with potential approximation errors. As the index increases, the divergence in results between these two methods becomes more significant, highlighting the importance of understanding the underlying mathematical principles.

Hope this helps.

[raven]

Thank you for your answer, I understand what you mean, it's just that I have to use the expansion method in my research, but this error is too large, so I want to find some way to solve this problem

[Umar]

Hi @Raven,

No problem. I understand your concern. However, @Torsten and @Steven Lord provided some good information. Does any of that helps you out resolving your problem.

[raven]

I get the idea, thank you very much for your help