matlab练习2（多项式插值，数值积分，数据拟合的最小二乘法，常微分方程数值解）

不吃橘子的橘猫

已于 2023-05-13 12:31:38 修改

阅读量655

点赞数 2

分类专栏：数学建模文章标签： matlab python 开发语言

于 2023-05-13 02:00:00 首次发布

本文链接：https://blog.csdn.net/weixin_73011353/article/details/130589026

版权

数学建模专栏收录该内容

23 篇文章 1 订阅

订阅专栏

文章展示了多项式插值用于求解sin(π/8)的近似值及截断误差估计；接着应用复化梯形公式和复化Simpson公式进行数值积分；然后通过最小二乘法对化学反应中质量浓度与时间的数据进行拟合；最后，采用预估校正Euler法求解常微分方程的数值解并绘制图形。

摘要由CSDN通过智能技术生成

7.多项式插值: 由函数y=sin x在三点0，π/4，π/2处的函数值，构造二次插值多项式P2(x)，计算sin(π/8)的近似值，并估计截断误差。

x = [0, pi/4, pi/2];
y = sin(x);
f = (y(3) - y(2)) / (x(3) - x(2)) - (y(2) - y(1)) / (x(2) - x(1));
P2 = @(x) (3 - sqrt(2))/8 * sin(x) + (sqrt(2) - 1)/8 * sin(pi/4) + (3 + sqrt(2))/8 * sin(pi/2);
approximate_value = P2(pi/8);
truncation_error = f * (pi/8 - 0) * (pi/8 - pi/4) * (pi/8 - pi/2);
fprintf('sin(pi/8)的近似值为：%.4f\n', approximate_value);
fprintf('截断误差的估计值为：%.4f\n', truncation_error);

8.数值积分: 已知函数 $f(x)=\frac{sin(x)}{x}$ 的如下数据表：

1/8

1/4

3/8

1/2

5/8

3/4

7/8

f(x)

0.997

3978

0.989

6158

0.976

7267

0.958

8110

0.936

1556

0.908

8517

0.877

1926

0.841

4710

试分别用复化梯形公式和复化Simpson公式计算定积分 $I=\int_{0}^{1}\frac{sin(x)}{x}dx$ 的近似值

复化梯形公式:

a = 0; 
b = 1; 
n = 100; 
h = (b - a)/n; 
x = a:h:b; 
f = zeros(1,n+1); % 初始化
for i = 1:n+1 % 计算f
    if x(i) == 0 % 处理x等于0的情况
        f(i) = 1;
    else
        f(i) = sin(x(i))/x(i);
    end
end
s = sum(f) - (f(1) + f(n+1))/2; 
I = h*s;
disp(['使用复化梯形公式计算得到I的值为：', num2str(I)]);

复化Simpson公式:

a = 0;
b = 1;
n = 100;
dx = (b-a)/n;
x = a+dx:dx:b-dx;
y = sin(x)./x;
S = dx/3 * (sin(a)/(a+eps) + sin(b)/(b-eps) + 4*sum(y(1:2:end-1)) + 2*sum(y(2:2:end-2)));
disp(['使用复化Simpson公式得到I的值为：', num2str(S)]);

9.数据拟合的最小二乘法: 在某化学反应中，生成物的质量浓度y（单位：10-3g/cm3）与时间t（单位：min）的关系为 $y=\frac{t}{\alpha +\beta t}$ ，现测得一组数据如下：

t1	1	2	3	4	6	8	10	12	14	16
y1	4.00	6.41	8.01	8.79	9.53	9.86	10.33	10.42	10.53	10.61

t = [1 2 3 4 6 8 10 12 14 16];
y = [4.00 6.41 8.01 8.79 9.53 9.86 10.33 10.42 10.53 10.61];
z = 1 ./ y;
p = polyfit(t, z, 1);
b_over_a = p(1);
a = 1 / p(2);
b = b_over_a * a;
disp(['a = ', num2str(a)]);
disp(['b = ', num2str(b)]);
tt = linspace(0, 16, 100);
yy = a ./ (1 + b * tt);
plot(t, y, 'o', tt, yy);
xlabel('时间t(min)');
ylabel('质量浓度y(10^-3 g/cm^3)');
legend({'原始数据', '拟合曲线'});

试用数据拟合的最小二乘法确定上述经验公式中的参数 $\alpha$ 和 $\beta$ ，并绘制拟合曲线图。

10.常微分方程数值解: 用预估校正Euler法，求出步长h=0.1的所有点的值，并绘制图形。

$\left\{\begin{matrix}y{}'=y-\frac{2x}{y},x\epsilon [0,1] \\y(0)=1 \end{matrix}\right.$

f = @(x, y) y - 2*x/y;

x0 = 0;
y0 = 1;
h = 0.1;
xn = 1;
x = x0:h:xn;

for i = 1:length(x)-1
    y1 = y0 + h*f(x(i), y0);
    y2 = y0 + h/2*(f(x(i), y0) + f(x(i+1), y1));
    y0 = y2;
    y(i+1) = y0;
end

plot(x, y)
title('Solution with h=0.1')
xlabel('x')
ylabel('y')