MATLAB运用——多项式插值

最新推荐文章于 2024-03-14 22:33:18 发布

闲谈社

最新推荐文章于 2024-03-14 22:33:18 发布

阅读量4.3k

点赞数 2

分类专栏： MATLAB语言文章标签： matlab

本文链接：https://blog.csdn.net/Ricardo2/article/details/116790924

版权

MATLAB语言专栏收录该内容

23 篇文章 20 订阅

订阅专栏

实验研究了误差传播在积分序列中的影响，发现算法一随着迭代次数增加误差增大，而算法二表现更稳定。此外，通过拉格朗日插值与三次样条插值对比，揭示了高阶多项式插值可能导致Runge震荡，而样条插值能更好地逼近原函数。

摘要由CSDN通过智能技术生成

实验1-误差传播与算法稳定性

问题提出

考虑一个简单的由积分定义的序列
$\int_{0}^{1}x^ne^{x-1}dx，n=1,2,……$ ————式（1）
显然 $E_n>0(n=1,2,……)$ ，当n=1时，利用分部积分得
$E_1=1/e$
而对于n≥2，经分部积分可得
$E_n=\int_{0}^{1}x^ne^{x-1}dx=x^ne^{x-1}|_{0}^{1}-\int_{0}^{1}nx^{n-1}e^{n-1}ds，n=1,2,……$
也就是
$E_n=1-nE_{n-1}，n=2,3,……$ ————式（2）
而再由式（1），可得
$E_n=\int_{0}^{1}x^ne^{x-1}dx≤\int_{0}^{1}x^ndx=1/(n+1)$ ————式（3）

实验内容

由递推关系式（2）可以得到积分序列 ${E_n\}$ 的两种算法。
算法一：式（2）的直接应用，即
$E_1=1/e，E_n=1-nE_{n-1}，n=1,2,……$
算法二：利用式（2）变形得到 $E_N=0$
$E_{n-1}={1-E_n \over n}，n=N-1,N-2,……,3,2$

实验要求

从程序观察两种算法的稳定性，再从理论验证。

程序

function f2
%编写算法1和算法2的运行程序，得出序列值和误差
%绘图比较两种算法
Step=str2num(char(inputdlg('输入迭代步数n，要求n大于等于1','步数')));
if Step<1
    errordlg('输入迭代步数n，要求n大于等于1');
    return
end
Sd=str2num(char(inputdlg('输入有效数字','有效数字')));
format long;

%用库函数integral计算
for n=1:Step
    fun=@(x)x.^n.*exp(x-1);
    func(n)=integral(fun,0,1);
end

%算法一
digits(Sd);
func1(1)=vpa(1/exp(1));
for n=2:Step
    func1(n)=vpa(1-n*func1(n-1));
end
err1=abs(func1-func);

%算法二
digits(Sd);
func2(Step)=0;
for n=Step:-1:2
    func2(n-1)=vpa((1-func2(n)/n));
end
err2=abs(func2-func);

clf
subplot(1,2,1)
plot([1:Step],func,'b-');
hold on
plot([1:Step],func1,'g-');
grid on
hold on
plot([1:Step],func2,'r--');
xlabel('n');ylabel('En');
text(2,func1(2),'func1');text(4,func2(4),'func2');text(1,func2(1),'func');

subplot(1,2,2)
plot([1:Step],err1,'-');
grid on
hold on
plot([1:Step],err2,'r--');
text(2,err1(2),'err1');text(4,err2(4),'err2');
xlabel('n');ylabel('Err');

运行程序

当迭代步数为10，有效数字为7时
在这里插入图片描述

当迭代步数为50，有效数字为7时
在这里插入图片描述
可见，迭代次数比较小时，算法一的误差波动比较大。当迭代次数比较大时，算法一的误差会越大。相比下，算法二更为稳定。

理论

设算法一的 $E_1$ 的误差为 $e_1$ ，则 $E_n$ 的误差为
$e(E_n)=e(1-nE_{n-1})=ne(E_{n-1})=……=n^{n-1}e(E_1)=n^{n-1}e_1$
设算法一的 $E_N$ 的误差为 $ε _N$ ，由 $E_N$ 向前递推得到的 $E_n$ 的误差为 $ε _n$
$(E_{n-1})=ε({1-E_n \over n})={ε(E_n)\over n}={ε(n)\over n}，$

明显可见，当n越大，算法一的误差越大，而相反算法二的误差越小。

实验2-多项式插值的震荡和样条插值的收敛

问题提出

考虑在一个固定区间上用插值逼近一个函数，显然Lagrange插值中使用的节点越多，插值多项式的次数越高。那么，当插值多项式的次数增加时，插值多项式是否也更加逼近函数呢？下面给出一个例子模拟Runge震荡现象。同时也对比样条插值的收敛情况。

实验内容

设区间[-5,5]上函数
$f(x)={1\over{1+25x^2}}$
在区间上做等距划分，分点为
$x_i=-1+{2i\over n}，i=0,1,2,……,n$
则拉格朗日多项式为
$L_n(x)=\sum_{i=o}^{n}{1\over 1+25x_i^2}l_i(x)$
其中， $l_i(x),i=1,2,……,n$ ，是n次Lagrange插值基函数

实验要求

（1）选择不断增大的分点数数目，画出原函数 $f (x)$ ，插值多项式函数 $L_n(x)$ 及3次样条插值函数在[-5,5]上的图像
（2）选择其他函数，例如定义在区间[-5,5]上的函数
$h(x)={x\over 1+x^2}，g(x)=arctanx$
（3）样条插值最早产生于工业部门，考虑如下问题，某汽车制造商用3次样条插值设计车门的曲线，其中一段的数据如下所示
在这里插入图片描述

程序

（1）

function f2
%选择不同的函数和多项式次数，输出原函数、拉格朗日插值函数和三次样条插值函数的图像。

n=str2num(char(inputdlg('输入插值多项式的次数N:')));
Nb=char(inputdlg('请选择函数，f(x)、g(x)、h(x)分别对应f、g、h'));
switch Nb
    case 'f'
        func=@(x)1./(1+25*x.^2);
    case 'h'
        func=@(x)x./(1+x.^4);
    case 'g'
        func=@(x)atan(x);
end      

if n<1
    errordlg('插值多项式的次数输入错误');
    return;
end

x0=linspace(-5,5,n+1);
y0=feval(func,x0);
x=-5:0.1:5;
%拉格朗日插值
Lg=lagrange(x0,y0,x);
%样条插值插值
Sp=interp1(x0,y0,x,'spline');

%作图
clf;
subplot(1,2,1);
fplot(func,[-5,5],'r-');
hold on;
plot(x,Lg,'b--');
xlabel('x');ylabel('y');title('拉格朗日插值');
subplot(1,2,2);
fplot(func,[-5,5],'r-');
hold on;
plot(x,Sp,'b--');
xlabel('x');ylabel('y');title('样条插值');


function Lg=lagrange(x0,y0,x)
%x0为节点的x坐标
%y0为节点的y坐标
%x为插值点的x坐标

if (nargin==3)%当输入参数为3个时，输出对应的插值
    n=length(x0);
    m=length(x);
    for i=1:m
        l=0;
        for k=1:n
            s=1;
            for j=1:n
                if k~=j
                    s=s*(x(i)-x0(j))/(x0(k)-x0(j));
                end
            end
            l=l+s*y0(k);
        end
        Lg(i)=l;
    end
elseif (nargin==2)%当输入参数为2个时，输出插值多项式
    syms x;
    n=length(x0);
    l=0;
    for k=1:n
        s=1;
        for j=1:n
            if k~=j
                s=s*(x-x0(j))/(x0(k)-x0(j));
            end
        end
        Lg=l+s*y0(k);
        Lg=collect(Lg);
        Lg = vpa(Lg,6);
    end
end

（2）

x0=0:10;
y0=[0.8,0.0,0.79,1.53,2.19,2.71,3.03,3.27,2.89,3.06,3.19,3.29,0.2];
x=0:0.1:10;
cs=spline(x0,y0);%当端点斜率已知时，可以用两个额外元素指定值向量 y（一个元素在起点，一个元素在终点）以定义端点斜率。
y=ppval(cs,x);%ppval计算在插值区间中各点上的样条拟合。
plot(x,y)