插值多项式的值与原函数值的最大偏差模数

Question

我正在学习 GNU Octave 的计算，并且正在努力完成以下任务。

给定插值节点X，定义函数lagrange(X,Y,x0)和original_function(x0)，分别计算插值多项式和原插值函数在x0点的值。并且已经定义了 delta_max(Y1,Y2) 函数。

我需要编写代码，查看区间 [xmin,xmax] 中的每个点，步长为 h=xmax−xmin / 100N（N 是插值区间的数量），并找到偏差的最大模数达到了从原始函数值插值多项式的值，以及这个最大值本身的值。

任务的输入/输出格式：

示例输入：

1 2 3

样本输出：

h = 0.010000
x_delta_max = 2.6100
delta = 0.58020

我很难在给定的时间内（13秒）提供答案。 有什么明显的我想念的吗？

我需要使用 GNU Octave 完成这项任务，但如果对代码的整体结构有任何建议（例如，我犯的错误不是特定于八度的），MatLab 示例将很有用。

首先，我计算了这个区间的步长和所有点的数组。 x_h 用于 X 坐标，y_h 用于 Y。 在下一步中，我尝试计算每个给定点中的插值多项式的值（如任务中所述）。我认为此时此刻我做错了什么，因为最后一个操作（delta_max）相对简单。

h = (X(end) - X(1)) / (100 * (size(X, 2) - 1));
x_h = X(1):h:X(end);
y_h = arrayfun(@original_function, x_h);
lagr = arrayfun(@(x) lagrange(x_h, y_h, x), x_h);
[delta, x_delta_max] = delta_max(y_h, lagr);

但是这段代码执行时间太长，而且我总是没有时间限制。

下面是 lagrange(X, Y, x0), original_function(x0) 和 delta_max(Y1, Y2) 的代码。

lagrange() 函数计算给定点的插值多项式的值。 X - 插值节点的x坐标，Y - 插值节点的y坐标，x0 - 计算值的点

function y0 = lagrange(X,Y,x0)
    y0 = 0;
    n = size(X, 2);
    for i = 1:n
        p = 1;
        for j = 1:n
            if (j != i)
                p *= (x0 - X(j)) / (X(i) - X(j));
            endif;
        endfor;
        y0 += Y(i) * p;
    endfor;
    return;
endfunction

original_function(x0)

function y0=original_function(x0)
    y0 = exp(x0);
endfunction

delta_max(Y1, Y2)

function [delta, i] = delta_max(Y1,Y2)
    delta = 0;
    i = 1;
    n = size(Y1, 2);
    for j=1:n
        if abs(Y1(j) - Y2(j)) > delta
            delta = abs(Y1(j) - Y2(j));
            i = j;
        endif;
    endfor;
endfunction

Answer 1

我做了一些修改和评论

clear

function y0 = lagrange(X,Y,x0)
    y0 = 0;
    n = size(X, 2);
    for i = 1:n
        p = 1;
        for j = 1:n
            if (j != i)
##            Use .* for elementwise operating
              p = p.*(x0 - X(j)) / (X(i) - X(j));
            endif;
        endfor;
        y0 += Y(i) * p;
    endfor;
## I removed the "return" command
endfunction

function y0=original_function(x0)
    y0 = exp(x0);
endfunction

function [delta, i] = delta_max(Y1,Y2)
    delta = 0;
    i = 1;
    n = size(Y1, 2);
    for j=1:n
        if abs(Y1(j) - Y2(j)) > delta
            delta = abs(Y1(j) - Y2(j));
            i = j;
        endif;
    endfor;
endfunction

# You can call the functions directly by its name
X = [1, 2, 3]
Y = original_function(X)
h = (X(end) - X(1)) / (100 * (size(X, 2) - 1))
x_h = X(1):h:X(end)
y_h = original_function(x_h)
lagr = lagrange(X, Y, x_h)
[delta, x_delta_max] = delta_max(y_h, lagr);