matlab中龙贝格积分,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

原文：https://www.cnblogs.com/AIxiaodi/p/matlab_0001_romberg.html