Matlab
% ----------------------------------------------------------------------------------------------- % Written by QQ742971636 % ----------------------------------------------------------------------------------------------- % 习题 对实验数据进行 lagrange 插值、三次样条插值 % 绘制出实验数据散点图,lagrange插值函数图,三次样条插值图 % ------------------------------------------------------------------------------------------------ % -------------------------------------------------------------------------------------------- % 实验数据 x y 初始赋值,画出原始数据折线图 % -------------------------------------------------------------------------------------------- clc,clear x=[0.0 0.5 1.0 6.0 7.0 9.0]; y=[0.0 1.6 2.0 1.5 1.5 0.0]; scatter(x,y,'k'); %画出原始数据散点图 hold on % -------------------------------------------------------------------------------------------- % lagrange插值过程 % -------------------------------------------------------------------------------------------- syms t; %系统变量 n = length(x); %x长度 f = 0.0; %lagrange插值多项式初始化 for i = 1:n %循环计算基函数 l = y(i); for j = 1:i-1 %计算基函数某个项前面一截,i之前的 l = l*(t-x(j))/(x(i)-x(j)); end for j = i+1:n %计算基函数某个项后面一截,i之后的 l = l*(t-x(j))/(x(i)-x(j)); end f = f + l; %得出最后的lagrange插值多项式 simplify(f); %化简多项式 if(i==n) f = collect(f); %合并同类项 f = vpa(f,6); %截断表达 end end disp("lagrange插值多项式结果:"); disp(f); %打印最后的lagrange插值多项式 % -------------------------------------------------------------------------------------------- % 画图 % -------------------------------------------------------------------------------------------- x0=linspace(0,9,200); f0 = subs(f,'t',x0); plot(x0,f0,'r'); %画图 hold on % -------------------------------------------------------------------------------------------- % 实验数据 x y 初始赋值 % -------------------------------------------------------------------------------------------- x=[0.0 0.5 1.0 6.0 7.0 9.0]; y=[0.0 1.6 2.0 1.5 1.5 0.0]; % -------------------------------------------------------------------------------------------- % 三次样条 插值过程 % -------------------------------------------------------------------------------------------- n = length(x); %计算长度 a = y(1 : end - 1); %插入n-1个 b = zeros(n - 1, 1); % 三次样条插值函数的系数 d = zeros(n - 1, 1); dx = diff(x); %计算剖分差距h dy = diff(y); %计算剖分差距 A = zeros(n); B = zeros(n, 1); A(1, 1) = 1; %边界条件 A(n, n) = 1; %边界条件 % -------------------------------------------------------------------------------------------- % 构建系数矩阵A*c=B % -------------------------------------------------------------------------------------------- for i = 2 : n - 1 A(i, i - 1) = dx(i - 1); A(i, i) = 2*(dx(i - 1) + dx(i)); A(i, i + 1) = dx(i); B(i) = 3*(dy(i) / dx(i) - dy(i - 1) / dx(i - 1)); end c = A \ B; %计算c % -------------------------------------------------------------------------------------------- % 三次样条插值函数的系数 % -------------------------------------------------------------------------------------------- for i = 1 : n - 1 d(i) = (c(i + 1) - c(i)) / (3 * dx(i)); b(i) = dy(i) / dx(i) - dx(i)*(2*c(i) + c(i + 1)) / 3; end % -------------------------------------------------------------------------------------------- % 将点带入插值公式计算结果 % -------------------------------------------------------------------------------------------- xx=linspace(0,9,200); XISHU=zeros(1,4); % 三次样条系数 [mm, nn] = size(xx); yy = zeros(mm, nn); for i = 1 : mm*nn % 计算每一个元素 for ii = 1 : n - 1 % 选取分段函数系数 if xx(i) >= x(ii) && xx(i) < x(ii + 1) j = ii; break; elseif xx(i) == x(n) j = n - 1; end end yy(i) = a(j) + b(j)*(xx(i) - x(j)) + c(j)*(xx(i) - x(j))^2 + d(j)*(xx(i) - x(j))^3;%计算插值结果 XISHU=[XISHU;[a(j) , b(j),c(j),d(j)]]; end XISHU=unique(XISHU,'rows'); % 分段函数系数 XISHU=XISHU(2:6,:); disp("三次样条分段函数系数矩阵:"); disp(XISHU); % -------------------------------------------------------------------------------------------- % 画三次样条图 % -------------------------------------------------------------------------------------------- plot(xx,yy,'b'); hold on legend("原始数据","lagrange插值","三次样条插值"); title("插值结果图")%标题
结果:
本文作者:Dong
本文链接:
版权声明:本博客所有文章除特别声明外,均采用 CC BY-NC。本作品采用《知识共享署名-非商业性使用 4.0 国际许可协议》进行许可。您可以在非商业用途下自由转载和修改,但必须注明出处并提供原作者链接。 许可协议。转载请注明出处!