多项式拟合之Arnoldi算法[个人记录用，不建议观看]

多项式拟合之Arnoldi算法

之前写过一篇多项式拟合的方程法，但是之前版本，包括matlab自带的多项式拟合函数，效果都不是特别好。在这里介绍一种多项式拟合法，由于本人能力有限，所以原理的关键地方也不懂。

多项式拟合一般表达式

$$\begin{array}{rl}y& =a_0+a_{1}x+a_2{x}^2+\cdots +a_n{x}^n\\ & =\sum _{i=0}^n{a}_i{x}^i\end{array}$$

构建方程组
$$\begin{array}{rl}{y}_1& =a_0+a_1{x}_1+a_2{x}_1^2+\cdots +a_n{x}_1^n\\ y_2& =a_0+a_1{x}_2+a_2{x}_2^2+\cdots +a_n{x}_2^n\\ ⋮& \\ y_m& =a_0+a_1{x}_m+a_2{x}_m^2+\cdots +a_n{x}_m^n\end{array}$$
令：
$$X=\left[\begin{array}{cccc}{x}_1& 0& \cdots & 0\\ 0& x_2& \cdots & 0\\ ⋮& ⋮& \cdots & ⋮\\ 0& 0& \cdots & x_m\end{array}\right]$$

$$b=\left(\begin{array}{c}{b}_0\\ b_1\\ ⋮\\ b_n\end{array}\right)$$

$$\begin{gathered} q_0=[1,\cdots ,1{]}^T \\ \begin{array}{rl}A& =[q_0,X{q}_0,X^2{q}_0,\cdots ,X^n{q}_0]\\ & =\left[\begin{array}{cccc}1& x_1& \cdots & x_1^n\\ 1& x_2& \cdots & x_2^n\\ ⋮& ⋮& \cdots & ⋮\\ 1& x_m& \cdots & x_m^n\end{array}\right]\end{array} \end{gathered}$$

$$Ab=f$$

Krylov 子空间

令$v\in R^N$。“,由$v,Av,\cdots ,A^{m-1}v$ 所生成的子空间称之为由$v$与$A$所生成的第$m$个Krylov子空间,并记$K_m(v,A)$。

剩余的推导

略，能力有限，见谅

代码

polyfitA.m

function [d,H] = polyfitA(x,f,n)
m = length(x);
Q = ones(m,1);
H = zeros(n+1,n);
for k = 1:n
q = x.*Q(:,k);
for j = 1:k
H(j,k) = Q(:,j)'*q/m;
q = q - H(j,k)*Q(:,j);
end
H(k+1,k) = norm(q)/sqrt(m);
Q = [Q q/H(k+1,k)];
end
d = Q\f;

polyvalA.m

function y = polyvalA(d,H,s)
M = length(s);
W = ones(M,1);
n = size(H,2);
for k = 1:n
w = s.*W(:,k);
for j = 1:k
w = w - H(j,k)*W(:,j);
end
W = [W w/H(k+1,k)];
end
y = W*d;

测试

测试代码

clear
close all
%% 准备数据
x = 1:0.1:50;
f=0.001*x.^3.2+0.04*x.^2+20*sin(x)-30*cos(x)+20*rand(size(x));
n=50;
%% matlab自带多项式拟合
c= polyfit(x',f',n);
%% 本文多项式拟合函数
[d,h]= polyfitA(x',f',n);

%% 测试效果
y1=polyval(c,x');% matlab自带
y2 = polyvalA(d,h,x'); %本文函数
err1 = sqrt(f*y1/length(f));
err2=sqrt(f*y2/length(f));
figure(1)
plot(x,f,x,y1','k') %matlab
title("matlab自带函数误差："+err1);
figure(2)
plot(x,f,x,y2','r') %本文
title("本文函数误差："+err2);

而且，当多项式次数比较高的时候，比如100次，matlab自带函数会出现错误！

这是80次的情况