西电数值分析实验作业matlab 第二章

西电数值分析实验作业matlab 第二章

gihub代码地址

第一题

拉格朗日

%% e21Lagrange.m
clear all
%% 数据输入
format long
x = [0.4 0.5 0.6 0.7 0.8];
y = [-0.916291 -0.693147 -0.510826 -0.356675 -0.223144];
n = length(x);

%% 插值
x0 = 0.65;

%% 插值计算
res = 0;
for i=1:n
   l = 1;
   for j = 1:n %计算lj
       if j ~= i
           l = l*(x0-x(j))/(x(i)-x(j));
       end
   end
   res = res + l*y(i);
end
res

%% 误差计算
%ln(x)的n阶导除以n的阶乘再乘(x0-x(i))连乘

%4阶误差估计
r1 = 1 / (5 * 0.4^5);%为了取最大可能误差，ξ取0.4
r2 = 1;
for i = 1:n
   r2 =r2 * abs(x0 - x(i));
end
R = r1*r2;
R %最大可能误差

%% 图形化输出
plot(x,y);
hold on
plot(x0,res, '*'); %插值

牛顿

%% e21Newton.m	

clear all
%% 数据输入
format long
x = [0.4 0.5 0.6 0.7 0.8];
y = [-0.916291 -0.693147 -0.510826 -0.356675 -0.223144];

n = length(x);
f = zeros(n, n + 1);%差商表

%% 插值
x0 = 0.65;

%% 插值计算
%差商表
f(1:n, 1) = x;
f(1:n, 2) = y;
for j = 3 : n + 1
    for i = j - 1 : n
    f(i, j) = (f(i, j-1) - f(i-1, j-1))/(f(i, 1) - f(i-j+2, 1));
    end
end
f

%插值多项式
res = y(1);
for i = 2:n
    t = 1;
    for j = 1:i-1
    t=t*(x0-x(j));
    end
    res = res + f(i,i+1) *  t;
end
res

%% 误差计算
%ln(x)的n阶导除以n的阶乘再乘(x0-x(i))连乘

%4阶误差估计
r1 = 1 / (5 * 0.4^5);%为了取最大可能误差，ξ取0.4
r2 = 1;
for i = 1:n
    r2 =r2 * abs(x0 - x(i));
end
R = r1*r2;
R %最大可能误差

%% 图形化输出
plot(x,y);
hold on
plot(x0,res, '*');%插值点

第二题

%% e22.m		
clear all
%% 数据输入
t = [0.25 0.5 1 1.5 2 3 4 6 8];
c = [19.21 18.15 15.36 14.10 12.89 9.32 7.45 5.24 3.01];

format long
x = [0.25 1.5 4,8];%采用三次插值
y = [19.21 14.1 7.45 3.01];

n = length(x);

%% 插值计算
for k = 1 : n
    V = 1;
    for i = 1 : n
        if k ~= i
            V = conv(V,poly(x(i))) / (x(k) - x(i));
        end
    end
    L1(k, :) = V; 
    l(k, :) = poly2sym(V);
end
Ls = y * l;
fprintf('多项式为:\nP(x)=%s\n',Ls);
for k=1:n
    L(k) =  L1(1,k)*y(1) + L1(2,k)*y(2)+ L1(3,k)*y(3)+L1(4,k)*y(4);
end

%% 画图
plot(t,c)
x = 0.25:0.1:8 ;
Ly = x.^3.*L(1)+x.^2.*L(2)+x.^1.*L(3)+L(4);
hold on
plot(x,Ly)%近似表达式

第三题

%% e23.m
clear all
%% 数据输入
format long
x = 0:0.4:2;
y = x.^7 - 1.2.*x.^5 + 2.3.*x.^4 + 2.3.*x.^3 - 5.6.*x + 1.9;

n = length(x);

%% 插值计算

for k = 1 : n
    V = 1;
    for i = 1 : n
        if k ~= i
            V = conv(V,poly(x(i))) / (x(k) - x(i));
        end
    end
    L1(k, :) = V; 
    l(k, :) = poly2sym(V);
end
Ls = y * l;
fprintf('多项式为:\nP(x)=%s\n',Ls);
for k=1:n
    L(k) =  L1(1,k)*y(1) + L1(2,k)*y(2)+ L1(3,k)*y(3) + L1(4,k)*y(4) + L1(5,k)*y(5) + L1(6,k)*y(6);
end
%% 画图
x = 1:0.1:2;
y = x.^7 - 1.2.*x.^5 + 2.3.*x.^4 + 2.3.*x.^3 - 5.6.*x + 1.9;
Ly = x.^5.*L(1)+x.^4.*L(2)+x.^3.*L(3)+x.^2.*L(4)+x.^1.*L(5)+L(6);
plot(x,y)
hold on
plot(x,Ly)%近似表达式

第四题

%% e24.m

clear all
%% 数据输入
format long
x = [1,31,61];
yc = [5*60+51, 5*60+17 5*60+10];%分钟
yl = [19*60+4, 19*60+38 19*60+50];
y = yl-yc;

n = length(x);

%% 插值计算
for k = 1 : n
    V = 1;
    for i = 1 : n
        if k ~= i
            V = conv(V,poly(x(i))) / (x(k) - x(i));
        end
    end
    L1(k, :) = V; 
    l(k, :) = poly2sym(V);
end
Ls = y * l;
fprintf('多项式为:\nP(x)=%s\n',Ls);
for k=1:n
    L(k) =  L1(1,k)*y(1) + L1(2,k)*y(2)+ L1(3,k)*y(3);
end
%% 求导
for k= 1:n
    dL(k) = L(k) * (n-k);
end
dx0 = -dL(2)/dL(1);
dx0 = round(dx0); %取整
Lydx0 = dx0.^2.*L(1)+dx0.^1.*L(2)+L(3);
%% 画图
plot(x,y,'*')
x = 1:1:61 ;
Ly = x.^2.*L(1)+x.^1.*L(2)+L(3);
hold on
plot(x,Ly)
plot(dx0,Lydx0,'x')%日照最长天与日照时间
mm = 5; %%月份
if dx0>31
    mm = 6;
    dx0 = dx0 -31; 
end
fprintf('日照最长天为：\n%d月%d日,\n与日照时间:\n%d分钟\n',mm,dx0,Lydx0);

参考：lagrange插值法:求拉格朗日插值多项式matlab实现(内附代码及例题)