Problem in finding the correct Chebyshev interpolation line

5 次查看（过去 30 天）
Ignatius Tommy Pratama
评论： Mathieu NOE 2024-4-2
Ran in:
Hi, I have problem in finding the correct Chebyshev interpolation line. I have inputted the calculation procedure in accordance with the Chebyshev theorem. However the regression line is apparently does not agree with the trend of the provided data point. Please let me know I missed something or how to improve my code. Your help is highly appreciated.
%This program is for governing a polinomial equation using Chebyshev
%interpolation formula
clc; clear all; close all;
t=tic;
%% Data input
x_data = [1830; 1920; 1940; 1975; 2000; 2020]; %x values
y_data = [0.089; 3.52; 4.86; 17.05; 25.5; 35.01]; %y values
%% Chebyshev interpolating formula
% Normalize x-values
x_min = min(x_data);
x_max = max(x_data);
x_data_norm = 2 * (x_data - x_min) / (x_max - x_min) - 1;
% Number of data points
n = length(x_data);
% Compute Chebyshev nodes
x_nodes = cos((2*(1:n) - 1) * pi / (2*n));
% Compute Chebyshev coefficients
a = zeros(n, 1);
for i = 1:n
T = cos((i-1) * acos(x_nodes));
a(i) = (2/n) * sum(y_data .* T');
end
% Define the Chebyshev polynomial function
T = @(k, x) cos((k-1) * acos(x));
% Evaluate the interpolating polynomial
x_query = linspace(-1, 1, 100); % Query points
y_interp = zeros(size(x_query));
for i = 1:n
y_interp = y_interp + a(i) * T(i, x_query);
end
%% Plotting the data and show polynomial equation
% Plot the given data points and the interpolating polynomial
% Construct the interpolation equation
interpolation_eqn = '';
for i = 1:n
if i == 1
interpolation_eqn = [interpolation_eqn num2str(a(i))];
else
interpolation_eqn = [interpolation_eqn ' + ' num2str(a(i)) ' * T' num2str(i-1) '(x)'];
end
end
disp(['Interpolation Equation:']);
Interpolation Equation:
disp(['P(x) = ' interpolation_eqn]);
P(x) = 28.6763 + -17.4761 * T1(x) + 3.8073 * T2(x) + -0.17701 * T3(x) + -0.17183 * T4(x) + -1.7569 * T5(x)
% Construct the interpolation equation in terms of powers of x
expanded_eqn = '';
for i = 1:n
if i == 1
expanded_eqn = [expanded_eqn num2str(a(i))];
else
expanded_eqn = [expanded_eqn ' + ' num2str(a(i)) ' * (2 * cos(' num2str(i-1) '* acos(x)) - 1)^' num2str(i-1)];
end
end
disp(['Expanded Interpolation Equation:']);
Expanded Interpolation Equation:
disp(['P(x) = ' expanded_eqn]);
P(x) = 28.6763 + -17.4761 * (2 * cos(1* acos(x)) - 1)^1 + 3.8073 * (2 * cos(2* acos(x)) - 1)^2 + -0.17701 * (2 * cos(3* acos(x)) - 1)^3 + -0.17183 * (2 * cos(4* acos(x)) - 1)^4 + -1.7569 * (2 * cos(5* acos(x)) - 1)^5
% Plot the data
figure;
plot(x_data, y_data, 'ro', 'MarkerSize', 10); % plot data points
hold on;
plot(x_query*(x_max-x_min)/2+(x_max+x_min)/2, y_interp, 'b-', 'LineWidth', 2); % plot interpolating polynomial
xlabel('Year');
ylabel('CO2 Concentration');
title('Chebyshev Interpolation');
legend('Given Data', 'Interpolating Polynomial', 'Location', 'best');
grid on;
%% Estimate P(x) for specific x values
% Evaluate the interpolating polynomial at specific x-values
x_values = [1915, 2010, 2050];
y_values = zeros(size(x_values));
% Evaluate the interpolating polynomial at the given x-values
for i = 1:length(x_values)
x_norm = 2 * (x_values(i) - x_min) / (x_max - x_min) - 1;
for j = 1:n
y_values(i) = y_values(i) + a(j) * T(j, x_norm);
end
end
disp(['For x = 1915, y = ', num2str(y_values(1))]);
For x = 1915, y = 27.4652
disp(['For x = 2010, y = ', num2str(y_values(2))]);
For x = 2010, y = 16.5339
disp(['For x = 2050, y = ', num2str(y_values(3))]);
For x = 2050, y = -30.148
toc(t)
Elapsed time is 0.692066 seconds.

请先登录，再回答此问题。

采纳的回答

Mathieu NOE
Mathieu NOE 2024-3-27
Ran in:
hello
maybe this ?
%% Data input
x_data = [1830; 1920; 1940; 1975; 2000; 2020]; %x values
y_data = [0.089; 3.52; 4.86; 17.05; 25.5; 35.01]; %y values
% Chebyshev polynomial interpolation
x = linspace(x_data(1),x_data(end),25);
% Barycentric interpolation
m = numel(x_data);
ts = -1+1/m:2/m:1-1/m;
c = (-1).^(0:m-1).*sin(pi*(ts+1)/2);
numer = zeros(size(x));
denom = zeros(size(x));
exact = zeros(size(x));
for j = 1:m
xdiff = x-x_data(j);
temp = c(j)./xdiff;
numer = numer + temp*y_data(j);
denom = denom + temp;
exact(xdiff==0) = j;
end
y = numer./denom;
jj = find(exact);
y(jj) = y_data(exact(jj));
% plot
plot(x_data,y_data,'*b',x,y,'r');

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Polynomials 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by