数值分析的几道上机编程题

原理都在书上,就不写了，很简单。

T1 四次牛顿插值

clc;
clear all;
close all;

%% 4次牛顿插值
x1=[0.2 0.4 0.6 0.8 1.0];
y1=[0.98 0.92 0.81 0.64 0.38];
len=length(y1);
a=y1; % 表第一列是f(x)
% 求均差ak,表
for j=2:len  
    for k=len:-1:j
        % 注意下标起点是1
        a(k)=(a(k)-a(k-1))/(x1(k)-x1(k-j+1));
    end
end
syms x df d;
df(1)=1;    % ak*xxx
d(1)=a(1);  % d=ak*xxx
% 求牛顿插值多项式的各项,公式3.2
for i=2:len 
    df(i)=df(i-1)*(x-x1(i-1));
    d(i)=a(i)*df(i);
end
% 求和得到牛顿插值多项式Pn(x)
P4=vpa(collect(sum(d)),5)         %P4即为4次牛顿插值多项式，并保留小数点后五位
figure;
% 绘制插值曲线
ezplot(P4,[0.2,1.08]);hold on
% 给出指定点 (xi,yi)
t=0.2+0.08.*[0,1,11,10];
plot(t,double(subs(P4,x,t)),'k*'); 
title('四次牛顿插值');grid on
legend('插值多项式P4(x)','指定点');

T3 8次多项式插值

clc;
clear all;
close all;

%% 8次多项式插值
x1=[0 1 4 9 16 25 36 49 64];
y1=[0 1 2 3 4 5 6 7 8];
n=length(y1);
syms x w Ln;
w=1;
% 求Wn+1(x),公式2.10
for i=1:n
    w=w*(x-x1(i));
end
% 求导
dw=diff(w,x);
% 求n次插值多项式,公式2.11
Ln=0;
for k=1:n
    Ln=Ln+w*y1(k)/(double(subs(dw,x,x1(k)))*(x-x1(k)));
end
% 显示插值多项式函数
L8=vpa(collect(Ln),5)
% 作图
figure;
ezplot(L8,[0,64]);hold on
plot(x1,y1,'k*');grid on
title('8次拉格朗日插值');
legend('插值多项式L8(x)','实际数据点');

# T16

clc;
clear all;
close all;

%% 原始数据
t=[0 0.9 1.9 3.0 3.9 5.0];
s=[0 10 30 50 80 110];
w=ones(1,length(t));
%% 绘制散点图
figure;
plot(t,s,'*');hold on
xlabel('时间 t/s');
ylabel('距离 s/m');
%% 观察各点在一条直线附近,拟合线性函数-
% s=a+bt-参考例9-a0=1,a1=t
a00=sum(w);
a01=sum(w.*t);
a11=sum(w.*t.^2);
a0s=sum(w.*s);
a1s=sum(w.*t.*s);
a=[a00 a01;a01 a11]\[a0s;a1s]
%% 绘制拟合的运动方程
t=0:0.1:5;
s=a(1)+a(2)*t;
plot(t,s,'k-');grid on
legend('实际数据','拟合方程');
%% 输出运动方程
syms s t;
s=vpa(a(1)+a(2).*t,6) 

T17

clc;
clear all;
close all;

%% 原始数据
x=[19 25 31 38 44];
y=[19.0 32.3 49.0 73.3 97.8];
w=ones(1,length(x));
%% 绘制散点图
figure;
plot(x,y,'*');hold on
xlabel('x');ylabel('y');
%% 指定拟合
% y=a+bx^2-参考例9-a0=1,a1=x^2
a00=sum(w);
a01=sum(w.*x.^2);
a11=sum(w.*x.^4);
a0y=sum(w.*y);
a1y=sum(w.*x.^2.*y);
a=[a00 a01;a01 a11]\[a0y;a1y]
%% 绘制拟合的运动方程
x1=19:0.1:44;
y1=a(1)+a(2)*x1.^2;
plot(x1,y1,'k-');grid on
legend('实际数据','拟合方程');
%% 输出运动方程
syms y1 x1;
y1=vpa(a(1)+a(2).*x1.^2,6) 
%% 计算均方误差
d=sqrt(sum(double(subs(y1,x1,x))-y).^2)

T18

clc;
clear all;
close all;

%% 原始数据
x=0:5:55;
y=[0 1.27 2.16 2.86 3.44 3.87 4.15 4.37 4.51 4.58 4.62 4.64].*10^(-4);
w=ones(1,length(x));
%% 绘制散点图
figure;
plot(x,y,'*');hold on
xlabel('时间 t/s');ylabel('浓度 y');
%% 观察各点在一条直线附近,拟合指数函数
% y=a*e^(-b/x)-参考例10
% 两边同时取对数 - yb=lna-b/t=A+bx - a0=1,a1=x 
% 提出(0,0)
w=w(2:end);
x=x(2:end);
y=y(2:end);
xb=-1./x;
yb=log(y);
a00=sum(w);
a01=sum(w.*xb);
a11=sum(w.*xb.^2);
a0yb=sum(w.*yb);
a1yb=sum(w.*xb.*yb);
a=[a00 a01;a01 a11]\[a0yb;a1yb]
%% 换回去
b(1)=exp(a(1));
b(2)=a(2);
b
%% 绘制拟合的运动方程
x1=0.1:0.1:55;
y1=b(1)*exp(-b(2)./x1);
plot(x1,y1,'k-');grid on
legend('实际数据','拟合方程');
%% 输出运动方程
syms y1 x1;
y1=vpa(b(1)*exp(-b(2)./x1),6)