【数值分析作业】

Experiment2_1

% 数值实验题2.1

clc,clear;
x = linspace(0.4, 0.9, 6);
ln_x = [-0.916291, -0.693147, -0.510826,...
    -0.356675, -0.223144, -0.105361];
x = x(3:6);
q_x = 0.65;
ln_x = ln_x(3:6);
lag_fx = Lagrange(x, ln_x, q_x);
fprintf('largange插值法结果为：%.5f \n', lag_fx);
new_fx = Newton(x, ln_x, q_x);
fprintf('Newton插值法结果为：%.5f \n', new_fx);
eplison_abs = abs(lag_fx - log(q_x));
eplison_rel = eplison_abs / abs(log(q_x));
fprintf('绝对误差为：%f;相对误差为：%f', eplison_abs, eplison_rel);

Experiment2_2

% 数值实验题2.2

clc,clear;
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];
fx = Lagrange(t, C)
% fprintf('C与t的关系式为:%s\n',fx)

Experiment2_3

% 数值实验题2.3
clc,clear;
x_n = linspace(0, 2, 6);
syms t
f(t) = t^7 - 1.2*t^5 + 2.3*t^4 +  2.3*t^3 - 5.6*t + 1.9;
y_n = f(x_n);
p_x = Lagrange(x_n, y_n)
fplot(f, 'r--', LineWidth=2, Marker='*')
hold on
fplot(p_x, 'b:', LineWidth=2,Marker='o')
legend('原函数f(x)','插值函数')
title('原函数f(x)和插值函数图像')

Experiment2_4

% 数值实验题2.3

clc,clear;
x = [1, 31, 61];  % 日期
y = [793, 861, 880];  % 时长，分钟
ft = Lagrange(x, y);
x0 = solve(diff(ft)==0)  % 求导数为0的点
ft_func = matlabFunction(ft);
y0 = ft_func(floor(x0));
% fprintf('驻点为：%.4f \n', x0);
is_max = double(diff(ft, 2)) < 0;
if(is_max)
    fprintf('在第 %d 天日照最长，最大值为: %.4f分钟 \n', floor(x0), y0)
else
    disp('不是最大值')
end

需要自定义函数

1. Lagrange.m

function L = Lagrange(x, y, x0)
%Lagrange 此处显示有关此函数的摘要
%   此处显示详细说明
%   x为已知数据点的x坐标
%   y为已知数据点的y坐标
%   x0为插值点的x坐标
%   f为求得的拉格朗日多项式
%   f0为在x0出的插值
syms t;                     %申明变量，否则不可使用，此处t相当于书上的x
if(length(x)==length(y))    %对输入数据检查，长度是否一致
    n=length(x);            %申明定义样本长度
else 
    disp('x与y的维数不相等，请检查输入数据')   
    return;   
end
L = 0;
for i = 1:n                      %从1到n，与25行构成累加
    e = y(i);                    %先取y值，再通过两段for循环使其与 基函数 相乘
    for j = 1:i-1
      e = e*((t-x(j))/(x(i)-x(j)));
    end
    for j = i+1:n
      e = e*((t-x(j))/(x(i)-x(j)));   
    end
    
    L = L + e;                   %累加计算拉格朗日插值函数
%     L = simplify(L);             %化简函数，其实在这个函数里面可有可无，先提一下，详细见图
    if(i == n)
        if(nargin == 3)          %当输入数字为三个的时候，计算插值点函数值
            L = subs(L,'t',x0);  %计算插值点的函数
            %解释一下subs的意思， 在函数L中用‘x0’替换‘t’，计算赋给L，详细见图
        else                     %输入两个的时候计算输出函数式
            L = collect(L);      %合并同类项，次数相等的系数合并，详细见图
            L = vpa(L,4);        %设置精度，为有效数字位数
        end
    end 
end

2. Newton.m

function  N  = Newton( x,y,t )
%NEWTON 牛顿插值法
% x:X坐标向量
% y:Y坐标向量
% t:代求插值点
syms   p ;   %定义符号变量
N = y(1);    %表示初始化为f（x0）
dd = 0;
dxs = 1;
n = length(x);  % 注意，这里的n与书上式子不一样，程序中n表示有n个值，书上式子中表示有n+1个值  
% 构造牛顿插值方法
for i = 1:n-1
    for j = i+1:n   % 两次循环嵌套，可以成功生成f[x0,x1]到f[x.,...,xn]
        dd(j) = (y(j)-y(i))/(x(j)-x(i)); % 注意大循环的最后一行，
        % 从第2次起，此处的y数组已经更新为上一阶的差商 
    end
    temp1(i) = dd(i+1);         %计算对应本次循环(x-x0).....(x-xi)部分的f[x0,x1,.....]
    dxs = dxs*(p-x(i));         %除f[x,x1,...]之外的(x-x0).....(x-xn-1);
    N = N + temp1(i)*dxs;       %累加得f（x）
    y = dd;                     %用于更新y的数组，使y数组变为变为这一阶的差商
end
  simplify(N);
%以上为计算部分，下面是输出规则；
if(nargin == 2)                 %当输入的参数为两个的时候，输出函数式
    N = subs(N,'p','x');         
    N = collect(N);              %合并同类项，次数相等的系数合并
    N = vpa(N,4);                 %设置精度，为有效数字位数
else
    %读取要插值点的向量长度，可以直接对多点插值计算，
    %表示如果最后一个参数输入一个数组的话，可以得到对应数量的结果
    m = length(t);
    for i = 1:m
       temp(i) = subs(N,'p',t(i));
    end
    N = temp;
end