实验三:Runge现象的产生和克服

程序1：lagrange多项式插值

syms f x p dp lx L;

f=1/(1+25*x^2);

N=input('请输入插值节点数N=');

xx=-1:2/N:1;

p=1; L=0;

ff=zeros(1,length(xx));

for i=1:(N+1)

    x=xx(i);

    ff(i)=eval(f);

    syms x;

    p=p*(x-xx(i));

end

dp=diff(p);

for j=1:(N+1)

    x=xx(j);

    k=eval(dp);

    syms x;

    lx=p/((x-xx(j))*k);

    L=L+lx*ff(j);

end

aa=[-0.96:0.1:-0.06,0,0.06:0.1:0.96];

for i=1:length(aa)

    x=aa(i);

    S(i)=eval(L);

    fff(i)=eval(f);

end

e=0;

for i=1:length(aa)

    e=e+(S(i)-fff(i))^2;

end

e=sqrt(e/(20*21));

fprintf('插值偏差为e=%.6f\n',e)

ezplot(f,[-1,1])

hold on

ezplot(L,[-1,1])

hold on

plot(xx,ff,'*')

hold on

plot(aa,S,'o')

hold off

程序2：分段线性插值

syms f x p lx;

f=1/(1+25*x^2);

N=input('请输入插值节点数N=');

xx=-1:2/N:1;

p=1; L=0;

ff=zeros(1,length(xx));

for i=1:(N+1)

    x=xx(i);

    ff(i)=eval(f);

end

syms x

for i=1:N

    for j=1:(N+1)

        if j==i

            lx(i,j)=(x-xx(i+1))/(xx(i)-xx(i+1));

        else if j==i+1

                lx(i,j)=(x-xx(i))/(xx(i+1)-xx(i));

            else

                lx(i,j)=0;

            end

        end

    end

end

p=lx*ff';

aa=[-0.96:0.1:-0.06,0,0.06:0.1:0.96];

for i=1:length(aa)

    x=aa(i);

    for j=1:N+1

        if x<xx(j)

            break

        end

    end

    S(i)=eval(p(j-1));

    fff(i)=eval(f);

end

e=0;

for i=1:length(aa)

    e=e+(S(i)-fff(i))^2;

end

e=sqrt(e/(20*21));

fprintf('插值偏差为e=%.6f\n',e)

ezplot(f,[-1,1])

hold on

xxx=(-1:0.01:1);

for i=1:length(xxx)

    x=xxx(i);

    for j=1:N+1

        if x<xx(j)

            break

        end

    end

    SS(i)=eval(p(j-1));

end

plot(xxx,SS,'r')

hold on

plot(xx,ff,'*')

hold on

plot(aa,S,'o')

hold off

程序3：三转角插值法

syms f x df s s1 s2 s3 s4;

f=1/(1+25*x^2);

df=diff(f);

N=input('请输入插值节点数N=');

h=2/N;

xx=-1:2/N:1;

p=1; L=0;

ff=zeros(1,length(xx));

for i=1:(N+1)

    x=xx(i);

    ff(i)=eval(f);

    dff(i)=eval(df);

end

syms x

for i=1:N

    s1=(x-xx(i+1))^2*(h+2*(x-xx(i)))*ff(i)/h^3;

    s2=(x-xx(i))^3*(h+2*(xx(i+1)-x))*ff(i+1)/h^3;

    s3=(x-xx(i+1))^2*(x-xx(i))*dff(i)/h^2;

    s4=(x-xx(i))^2*(x-xx(i+1))*dff(i+1)/h^2;

    s(i)=s1+s2+s3+s4;

end

aa=[-0.96:0.1:-0.06,0,0.06:0.1:0.96];

for i=1:length(aa)

    x=aa(i);

    for j=1:N+1

        if x<xx(j)

            break

        end

    end

    S(i)=eval(s(j-1));

    fff(i)=eval(f);

end

e=0;

for i=1:length(aa)

    e=e+(S(i)-fff(i))^2;

end

e=sqrt(e/(20*21));

fprintf('插值偏差为e=%.6f\n',e)

ezplot(f,[-1,1])

hold on

xxx=(-1:0.01:1);

for i=1:length(xxx)

    x=xxx(i);

    for j=1:N+1

        if x<xx(j)

            break

        end

    end

    SS(i)=eval(s(j-1));

end

plot(xxx,SS,'r')

hold on

plot(xx,ff,'*')

hold on

plot(aa,S,'o')

hold off

程序4：三弯矩插值法：

syms f x ddf s;

f=1/(1+25*x^2);

ddf=diff(diff(f));

N=input('请输入插值节点数N=');

h=2/N;

xx=-1:2/N:1;

p=1; L=0;

ff=zeros(1,length(xx));

for i=1:(N+1)

    x=xx(i);

    ff(i)=eval(f);

    ddff(i)=eval(ddf);

end

syms x

for i=1:N

    A=(ff(i+1)-ff(i))/h-h*(ddff(i+1)-ddff(i))/6;

    B=ff(i)-h^2*ddff(i)/6;

    s(i)=(xx(i+1)-x)^3*ddff(i)/(6*h)+(x-xx(i))^3*ddff(i+1)/(6*h)+A*(x-xx(i))+B;

end

aa=[-0.96:0.1:-0.06,0,0.06:0.1:0.96];

for i=1:length(aa)

    x=aa(i);

    for j=1:N+1

        if x<xx(j)

            break

        end

    end

    S(i)=eval(s(j-1));

    fff(i)=eval(f);

end

e=0;

for i=1:length(aa)

    e=e+(S(i)-fff(i))^2;

end

e=sqrt(e/(20*21));

fprintf('插值偏差为e=%.6f\n',e)

ezplot(f,[-1,1])

hold on

xxx=(-1:0.01:1);

for i=1:length(xxx)

    x=xxx(i);

    for j=1:N+1

        if x<xx(j)

            break

        end

    end

    SS(i)=eval(s(j-1));

end

plot(xxx,SS,'r')

hold on

plot(xx,ff,'*')

hold on

plot(aa,S,'o')

hold off

图1：观察Runge现象

图2：分段线性插值

图3：三转角插值：

图4：三弯矩插值：