Problem in finding the correct Chebyshev interpolation line

2 次查看(过去 30 天)
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
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 CenterFile Exchange 中查找有关 Interpolation 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by