Matlab将输入的数学函数声明为内联(例：龙贝格积分)

编写数值分析实验的小程序时，其中一个实验需要输入被积函数，之后使用龙贝格积分法进行数值积分。

一开始使用简单的方法,每次计算函数值时使用，eval(subs(fun,x,1))这种方法，但是计算会的很慢。
百度之后发现可以直接把字符串声明为内联函数，如下:

    fun = input('输入被积函数(例如x+1): ', 's');
     % 将fun声明为内联函数
    fun = inline(fun);

之后输入积分区间，精度。利用公式不断计算数值定积分即可。

整个程序如下:

% 🐉贝格积分法
% 让用户输入函数及积分区间
% 并根据输入的精度结束循环
% 最多计算到T[9,9](在matlab里是T(10, 10));
clear all;
clc;

syms x;

while (1)
    % 定义max m = 9, max k = 9
    T = zeros(10, 10);
    fun = input('输入被积函数(例如x+1): ', 's');
    % 将fun声明为内联函数
    fun = inline(fun);
    s = input('输入积分区间(例如s = [0, 1]): ','s');
    e = input('输入精度e: ');
    eval(s);
    a = s(1);
    b = s(2);
    
    T(1, 1) = (b - a) / 2 * (fun(a) + fun(b));
    
    % 处理类似sinx/x的情况
    if (isnan(T(1, 1)))
        T(1, 1) = (b - a) / 2 * (1 + fun(b));
    end
    T = calc_fline(T, fun, a, b);   % 将内联函数传递给函数
    T = calc_remain(T, e);
    T
    
    flag = input('是否继续(0退出): ');
    if (~flag)
        break;
    end
    
end

% 计算romberg数表的第一行
function T = calc_fline(T, fun, a, b)
    for i = 2:1:10
        temp = 0;
        for j = 1:1:(2^(i-2))
            temp = temp + fun(...
                (a + (2 * j - 1) * (b - a) / 2 ^ (i - 1)));
        end
        T(1, i) = 1 / 2 * (T(1, i-1) + (b - a) / 2 ^ (i - 2) * temp);
    end
end

% 计算第2行至第10行的romberg积分
function T = calc_remain(T ,e)
    for m = 2:1:10
        for k = 1:1:(11-m)  % 写成m-1出错
            T(m, k) = (4 ^ (m - 1) * T(m-1, k+1) - T(m-1, k)) ...
                / (4 ^ (m - 1) - 1);
            abs(T(m, 1) - T(m - 1,1))
            if (abs(T(m, 1) - T(m - 1, 1)) < e)
                break;
            end
        end
        if (T(m, 2) == 0)
            break;
        end
    end
end