Matlab





```
% -----------------------------------------------------------------------------------------------
% Written by QQ742971636
% -----------------------------------------------------------------------------------------------
% f(x)=(10/x)^2*sin(10/x);
% 取epsilon=10^(-4) ，计算[1,3]上的romberg积分
% ------------------------------------------------------------------------------------------------
a = 1;
b = 3;
epsilon = 1e-4;
f = @(x)(10/x)^2*sin(10/x);  
[R,k,T] = aa(f,a,b,epsilon);%调用函数
fprintf("积分结果是    %f\r\n",R); %打印积分结果
fprintf("k=  %f\r\n",k); %打印迭代次数
disp(T); %打印T-表
function [R,k,T]=aa(fun,a,b,tol)
% 龙贝格(Romberg数值求解公式)
% inputs:
%   -fun：积分函数句柄
%   -a/b：积分上下限
%   -tol：积分误差
% Outputs:
%   -R：Romberg积分值
%   -k：迭代次数
%   -T：整个迭代过程
k=0; % 迭代次数
n=1; % 区间划分个数
h=b-a;
T=h/2*(fun(a)+fun(b)); %计算T(0,0)
err=1; %积分误差初始化
while err>=tol
    k=k+1;
    h=h/2;
    tmp=0;
    for i=1:n
        tmp=tmp+fun(a+(2*i-1)*h);
    end
    T(k+1,1)=T(k)/2+h*tmp; %外推算法
    for j=1:k
        T(k+1,j+1)=T(k+1,j)+(T(k+1,j)-T(k,j))/(4^j-1); %外推算法
    end
    n=n*2;
    err=abs(T(k+1,k+1)-T(k,k)); %更新积分误差
end
R=T(k+1,4);
end
```


结果

```
>> Untitled
积分结果是    -1.426025
k=  7.000000
  -56.5195         0         0         0         0         0         0         0
  -52.2329  -50.8040         0         0         0         0         0         0
  -23.8564  -14.3976  -11.9705         0         0         0         0         0
   -6.8278   -1.1516   -0.2685   -0.0828         0         0         0         0
   -2.6815   -1.2994   -1.3093   -1.3258   -1.3307         0         0         0
   -1.7327   -1.4164   -1.4242   -1.4260   -1.4264   -1.4265         0         0
   -1.5022   -1.4254   -1.4260   -1.4260   -1.4260   -1.4260   -1.4260         0
   -1.4450   -1.4260   -1.4260   -1.4260   -1.4260   -1.4260   -1.4260   -1.4260
```

matlab romberg积分 龙贝格积分

Matlab





```
clc,clear
x=[0,1,2,3,4];
%y=[0,0.075,0.086,0.052,0.0024];
y=[0,0.037,0.025,0.037,0];
scatter(x,y,'k');                                                      %画出原始数据散点图
hold on
duoxiangshi = polyfit(x,y,4); %3次拟合 降幂排列
disp(duoxiangshi);
xi = 0:0.001:4.0;
yi = polyval(duoxiangshi,xi);   %多项式求值
plot(xi,yi,'b');         %观测数据点
hold on
xlabel('x');
ylabel('y');
legend('原始数据','拟合曲线');
```

结果

![在这里插入图片描述](/static/img/images/02781.png)

MATLAB 对试验点进行三次拟合

Matlab





```
% -----------------------------------------------------------------------------------------------
% Written by QQ742971636
% -----------------------------------------------------------------------------------------------
% 习题  对实验数据进行 lagrange 插值、三次样条插值
%          绘制出实验数据散点图，lagrange插值函数图，三次样条插值图
% ------------------------------------------------------------------------------------------------
% --------------------------------------------------------------------------------------------
% 实验数据 x y 初始赋值，画出原始数据折线图
% --------------------------------------------------------------------------------------------
clc,clear
x=[0.0 0.5 1.0 6.0 7.0 9.0];
y=[0.0 1.6 2.0 1.5 1.5 0.0];
scatter(x,y,'k');                                                      %画出原始数据散点图
hold on
% --------------------------------------------------------------------------------------------
% lagrange插值过程
% --------------------------------------------------------------------------------------------
syms t;                                                                 %系统变量
n = length(x);                                                       %x长度
f = 0.0;                                                                %lagrange插值多项式初始化
for i = 1:n                                                            %循环计算基函数
    l = y(i);
    for j = 1:i-1                                                       %计算基函数某个项前面一截，i之前的
        l = l*(t-x(j))/(x(i)-x(j));
    end
    for j = i+1:n                                                      %计算基函数某个项后面一截，i之后的
        l = l*(t-x(j))/(x(i)-x(j));
    end
    f = f + l;                                                           %得出最后的lagrange插值多项式
    simplify(f);                                                       %化简多项式
    if(i==n)
        f = collect(f);                                                %合并同类项
        f = vpa(f,6);                                                  %截断表达
    end
end
disp("lagrange插值多项式结果：");  
disp(f);                                                                %打印最后的lagrange插值多项式
% --------------------------------------------------------------------------------------------
% 画图
% --------------------------------------------------------------------------------------------
x0=linspace(0,9,200);
f0 = subs(f,'t',x0);
plot(x0,f0,'r');                                                       %画图
hold on
% --------------------------------------------------------------------------------------------
% 实验数据 x y 初始赋值
% --------------------------------------------------------------------------------------------
x=[0.0 0.5 1.0 6.0 7.0 9.0];
y=[0.0 1.6 2.0 1.5 1.5 0.0];
% --------------------------------------------------------------------------------------------
% 三次样条 插值过程
% --------------------------------------------------------------------------------------------
n = length(x);                                                       %计算长度
a = y(1 : end - 1);                                                 %插入n-1个
b = zeros(n - 1, 1);                                               % 三次样条插值函数的系数
d = zeros(n - 1, 1);
dx = diff(x);                                                          %计算剖分差距h
dy = diff(y);                                                          %计算剖分差距
A = zeros(n);
B = zeros(n, 1);
A(1, 1) = 1;                                                          %边界条件
A(n, n) = 1;                                                          %边界条件
% --------------------------------------------------------------------------------------------
% 构建系数矩阵A*c=B
% --------------------------------------------------------------------------------------------
for i = 2 : n - 1
    A(i, i - 1) = dx(i - 1);
    A(i, i) = 2*(dx(i - 1) + dx(i));
    A(i, i + 1) = dx(i);
    B(i) = 3*(dy(i) / dx(i) - dy(i - 1) / dx(i - 1));
end
c = A \ B;                                                             %计算c
% --------------------------------------------------------------------------------------------
% 三次样条插值函数的系数
% --------------------------------------------------------------------------------------------
for i = 1 : n - 1
    d(i) = (c(i + 1) - c(i)) / (3 * dx(i));
    b(i) = dy(i) / dx(i) - dx(i)*(2*c(i) + c(i + 1)) / 3;
end
% --------------------------------------------------------------------------------------------
%  将点带入插值公式计算结果
% --------------------------------------------------------------------------------------------
xx=linspace(0,9,200);
XISHU=zeros(1,4);                                               %  三次样条系数
[mm, nn] = size(xx);
yy = zeros(mm, nn);
for i = 1 : mm*nn                                                  %  计算每一个元素
    for ii = 1 : n - 1                                                 %  选取分段函数系数
        if xx(i) >= x(ii) && xx(i) < x(ii + 1)
            j = ii;
            break;
        elseif xx(i) == x(n)
            j = n - 1;
        end
    end
    yy(i) = a(j) + b(j)*(xx(i) - x(j)) + c(j)*(xx(i) - x(j))^2 + d(j)*(xx(i) - x(j))^3;%计算插值结果
    XISHU=[XISHU;[a(j) , b(j),c(j),d(j)]];
end
XISHU=unique(XISHU,'rows');                              %  分段函数系数
XISHU=XISHU(2:6,:);
disp("三次样条分段函数系数矩阵：");
disp(XISHU);
% --------------------------------------------------------------------------------------------
%  画三次样条图
% --------------------------------------------------------------------------------------------
plot(xx,yy,'b');
hold on
legend("原始数据","lagrange插值","三次样条插值");
title("插值结果图")%标题
```


结果：

![在这里插入图片描述](/static/img/images/02782.png)

MATLAB lagrange插值 三次样条插值

Matlab





```
clear
% creat all input variable
disp('（1）中为A1   （2）中为A2');
A1 = [3.01,6.03,1.99;
        1.27,4.16,-1.23;
        0.987,-4.81,9.34];
A2 = A1;
A2(1,1) = 3.00;
A2(3,1) = 0.990;
b = [1,1,1]';
% display
disp('MAT A1 is ');
disp(A1);
disp('MAT A2 is ');
disp(A2);
disp('MAT b is ');
disp(b);
disp('detA1 is ');
disp(det(A1));
disp('detA2 is ');
disp(det(A2));
disp('cond(A1,1) is ');
disp(cond(A1,1));
disp('cond(A1,2) is ');
disp(cond(A1,2));
disp('cond(A1,inf) is ');
disp(cond(A1,inf));
%列主元消去法          
disp('在（1）中的情况：');
A = A1;
B = [A,b];
n = length(b);
rankA=rank(A);
rankB=rank(B);
if rankB-rankA>0
    disp('此方程组无解.\n')
end
if rankA==rankB
    if rankA==n
        fprintf('rankA=rankB=%d,有唯一解.\n',n)
        X=zeros(n,1);
        for p=1:n-1
              t=find(abs(B(p:end,p))==max(abs(B(p:end,p))))+p-1;
              if abs(B(t,p))~=abs(B(p,p))
                    l=B(t,:);
                    B(t,:)=B(p,:);
                    B(p,:)=l;
              end  %列主元判断
              for k=p+1:n     
                    m= B(k,p)/ B(p,p);
                    B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);
              end
        end
        %把方程组系数矩阵A化为同解的上三角矩阵
        b=B(1:n,n+1);
        A=B(1:n,1:n); 
        X(n)=b(n)/A(n,n);
        for q=n-1:-1:1
            X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);
        end
        %从xn至x1逐个求解上三角方程组
    else
        disp('rankA=rankB<n，有无穷多解.')
        return
    end
end
disp('方程组的解为：');
disp(X)
%列主元消去法          
disp('在（2）中的情况：');
A = A2;
B = [A,b];
n = length(b);
rankA=rank(A);
rankB=rank(B);
if rankB-rankA>0
    disp('此方程组无解.\n')
end
if rankA==rankB
    if rankA==n
        fprintf('rankA=rankB=%d,有唯一解.\n',n)
        X=zeros(n,1);
        for p=1:n-1
              t=find(abs(B(p:end,p))==max(abs(B(p:end,p))))+p-1;
              if abs(B(t,p))~=abs(B(p,p))
                    l=B(t,:);
                    B(t,:)=B(p,:);
                    B(p,:)=l;
              end  %列主元判断
              for k=p+1:n     
                    m= B(k,p)/ B(p,p);
                    B(k,p:n+1)= B(k,p:n+1)-m* B(p,p:n+1);
              end
        end
        %把方程组系数矩阵A化为同解的上三角矩阵
        b=B(1:n,n+1);
        A=B(1:n,1:n); 
        X(n)=b(n)/A(n,n);
        for q=n-1:-1:1
            X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);
        end
        %从xn至x1逐个求解上三角方程组
    else
        disp('rankA=rankB<n，有无穷多解.')
        return
    end
end
disp('方程组的解为：');
disp(X)
%use matlib 
%disp('使用matlib直接求解');
%disp(A2\b);
```



结果

```
>> Untitled
（1）中为A1   （2）中为A2
MAT A1 is 
    3.0100    6.0300    1.9900
    1.2700    4.1600   -1.2300
    0.9870   -4.8100    9.3400
MAT A2 is 
    3.0000    6.0300    1.9900
    1.2700    4.1600   -1.2300
    0.9900   -4.8100    9.3400
MAT b is 
     1
     1
     1
detA1 is 
   -0.0305
detA2 is 
   -0.4070
cond(A1,1) is 
   5.5228e+04
cond(A1,2) is 
   3.0697e+04
cond(A1,inf) is 
   5.6751e+04
在（1）中的情况：
rankA=rankB=3,有唯一解.
方程组的解为：
   1.0e+03 *
    1.5926
   -0.6319
   -0.4936
在（2）中的情况：
rankA=rankB=3,有唯一解.
方程组的解为：
   56.3415
  -22.1547
  -17.3024
```

MATLAB 列主元消去法 解方程组

Matlab





```
clc
clear
n=input('请输入矩阵的阶数n：');
A=zeros(n);
%构造矩阵A
for m=1:n
    A(m,m)=20;
end
for m=1:n-1
    A(m,m+1)=-8;
    A(m+1,m)=-8;
end
for m=1:n-2
    A(m,m+2)=1;
    A(m+2,m)=1;
end
disp(A)
%%%%%%%%%%%%%%Jacobi迭代法求解
%用矩阵X的行向量来存放每一次迭代的解
for m=1:n
    X(1,m)=1;
    Y(1,m)=0;
end
fs=1;  %设置初始值，以便开始下方的循环
t=1;
while fs >10^(-6)
    for p=1:n
        Y(t,:)=X(t,:);   %第t次迭代
        Y(t,p)=0;   %与解有关，与第几次迭代无关
        X(t+1,p)=-Y(t,:)*A(p,:)'/A(p,p);   %第t次迭代的第p个解
    end
    C=zeros(1,n);
    for u=1:n
        C(u)=abs(X(t+1,u));
    end
    fs=max(C);
    t=t+1;
end
t-1      %显示迭代次数
```

MATLAB Jacobi迭代法求解

我是一个专注技术分享的博主,尤其爱搞深度学习和单片机电子。你现在看到的是我的个人博客网站。