数值微分&数值积分Matlab程序分享

本文摘要（由AI生成）：

本文介绍了利用三次样条拟合实验数据，并通过微分计算得到反映速率的方法。文章还讨论了数值积分在无法获得解析解时的应用，包括等距节点的牛顿-柯特斯求积公式和不等距节点的高斯型求积公式。文章展示了使用trapz函数和quad函数进行数值积分的便捷性，并比较了不同方法的拟合效果。通过数值运算中的微分和积分运算，可以灵活处理数据并得出有用结论。

    过冷水在日常学习中，经常利用实验数值做微分，积分工作，数值微积具有较强的实用性。 微积分的数值方法,不同于高等数学中的解析方法,尤其适合求解没有或很难求出微分或积分表达式的实际问题计算。本期过冷水给大家较全面的介绍一下数值微分和数值积分求解理论和 MATLAB 求解数值微分和数值积分方法。

    数值微分和数值积分与插值和拟合往往是密不可分的。如在进行数值微分时,往往针对的是离散数据点,利用插值和拟合常可以减少误差,保证结果的精确性。而数值积分基本思路也来自于插值法,如所积函数的形式比较复杂,可通过构造一个插值多项式来代替原函数,从而使问题大大简化。

数值微分

    解析法中函数的导数是通过取极限来计算。当函数形式不明确,如以表格给出自变量和因变量关系时就不能通过解析法求导数,实际工作中通常是需要求列表函数在节点和非节点处的导数值,这正是数值微分所要解决的问题。数值微分方法可近似求出某点的导数值,或者将函数在某点的导数用该点附近节点上的函数近似表示。比如在分子动力学中想要锁定研究某个原子的速度。

常用以下三种思路建立数值微分公式。

（1）从微分定义出发,通过差分近似处理得到数值微分近似解;

（2）从插值近似公式出发,对插值公式求导得到数值微分近似解;

(3）先用最小二乘拟合方法根据已知数据得到近似函数,再对此近似函数求微分得到数值微分的近似解。

过冷水详细给大家讲讲根据这三种思路出发，给出具体操作方法，采用MATLAB函数解决数值微分计算方法。

差分近似微分

在微积分中,根据导数的定义可知,一阶微分的计算可以对两个相邻点x h和x之间函数取极限求得:

取其达到极限前的形式,就得到以下微分的差分近似式:

式中三种不同表示形式依次是一阶前向差分、一阶后向差分和一阶中心差分来近似表示微分,其中一阶中心差分的精度较高。高阶微分项可以利用低阶微分项来计算,如一个二阶微分可以表示为

所对应的差分式有:

当用差分近似微分时,通常需要步长h较短,因此实际差分近似为除以一个小的数的运算,计算误差将会被放大,在使用时需注意。

差分的MATLAB实现

在MATLAB中,可用diff函数求向量相邻元素的差,diff(y)./ diff(x)则表示一阶前向差分。diff函数的调用形式为

Y=diff(x)
Y=diff(x,n)
Y=diff(x,n,dim)

X表示求差变量,可以是向量或矩阵,如果X为向量,则diff(X)表示相邻元素的差，即

[x(2)-x(1),x(3)-x(2),...,x(n)-x(n-1)]

如果X为矩阵,则diff(X)表示各列相邻元素的差，即

[x(2,:)-x(1,:),x(3,:)-x(2,:),...,x(n,:)-x(n-1,:)]

n表示diff 函数循环运算n次;

(3) 当X是矩阵时,dim可以指定diff函数作用的维,当dim=2时,表示行元素进行diff 函数的运算。

差分法验证

clear
x=linspace(1,10,400);
y=sin(x);
dy=cos(x);
dyy=diff(y)./diff(x);
x1=linspace(1,10,40);
y1=sin(x1);
dy1=diff(y1)./diff(x1);
figure1 = figure;
axes1 = axes('Parent',figure1);
hold(axes1,'on');
plot(x,dy,'DisplayName','解析式');
plot(x(2:end),dyy,'DisplayName','变量数400','LineWidth',3,'LineStyle','--');
plot(x1(2:end),dy1,'DisplayName','变量数40','LineWidth',3);
ylabel('cos(x)');
xlabel('x');
title('sin(x)微分函数');
box(axes1,'on');
hold(axes1,'off');
set(axes1,'FontSize',15,'LineWidth',3);
legend(axes1,'show');

三次样条插值函数求微分

    若三次样条插值函数S(x)收敛于f(x),那么导数S'(x)收敛于f'(x),因此用样条插值函数S(x)作为f(x)近似函数,不但彼此函数值非常接近﹐而且导数值也很接近。用三次样条插值函数求数值导数是可靠的,这是求数值微分的有效方法。

MATLAB求离散数据的三次样条插值函数微分方法分以下三个步骤。

(1) 对离散数据用spline或pchip得到到其三次样条插值函数pp;

(2) 可用fnder函数求三次样条插值函数的导数,其调用形式为dp= fnder( pp，dorder)

(3）可用fnval 函数求导函数在未知点处的导数值,其调用形式为v= fnval( dp，x)

 利用三次样条插值求微分的方法分别计算t=2 min，3 min，4 min及t=1.5 min,2.5 min，4.5 min时的降温速平dT/dt。

三次样条插值函数求数值微分的程序如下。

t=0:5;
T=[92.0 85.3 79.5 74.5 70.2 67.0];
pp=spline(t,T)
plot(t,T,'o');
hold on
fnplt(pp)
dp=fnder(pp)
t1=[2,3,4,1.5,2.5,4.5];
dT=fnval(dp,t1)
disp('不同时刻的降温速率：')
disp([t1;dT])
pp = 
  包含以下字段的 struct:
      form: 'pp'
    breaks: [0 1 2 3 4 5]
     coefs: [5×4 double]
    pieces: 5
     order: 4
       dim: 1
dp = 
  包含以下字段的 struct:
      form: 'pp'
    breaks: [0 1 2 3 4 5]
     coefs: [5×3 double]
    pieces: 5
     order: 3
       dim: 1
dT =
   -5.3722   -4.6722   -3.8389   -5.7972   -4.9889   -3.2222
不同时刻的降温速率：
    2.0000    3.0000    4.0000    1.5000    2.5000    4.5000
   -5.3722   -4.6722   -3.8389   -5.7972   -4.9889   -3.2222

过冷水平常用使用插值做数据提供，导数值的使用较少，使用dT=fnval(dp,t1)的场合较少，可以灵活变化用于求插值p=fnval(pp,t1)

fnval(pp,t1)
ans =
   79.5000   74.5000   70.2000   82.2917   76.9125   68.4292

最小二乘法拟合函数求微分

    在实际工程实验中实验观测的离散数据不可避免地含有随机误差时,用插值公式spline求数值微分虽然样本点处误差较小,但可能会使非样本点处产生较大误差,不能较好的反映出整体趋势。为此,可用最小二乘法拟合实验数据,获得函数模型,然后再对其求导数。由于拟合不要求曲线经过全部的数据点,这样处理求导结果有助于消除实验误差带来的影响。

        多项式具有良好的计算性质,是最常选用用于拟合的函数形式。除了多项式拟合，过冷水在实践中也经常使用高斯多项式，或自己给出能够描述实验观测点趋势的函数。本案例过冷水给大家讲多项式拟合，利用MATLAB提供的多项式拟合函数polyfit和三次样条拟合函数csaps的拟合结果可以方便地求微分。

MATLAB利用多项式拟合函数polyfit求微分涉及以下三个步骤。

(1) 对离散数据用polyfit函数得到多项式系数向量p;

(2)用polyder函数求多项式拟合函数的导数;

(3)用polyval 函数求导函数在未知点处的导数值。

 系统的动力学模型为,根据表中数据确定模型参数m和k参数拟合需要反应速率,采用三次样条拟合求微分方法计算；系统的动力学模型为非线性形式，将其线性化，拟合后作为参数初始值，最后采用nlinfit函数直接拟合。

t=[0 20 40 60 120 180 300];
CA=[10 8 6 5 3 2 1];
plot(t,CA,'r*')
hold on
pp=csaps(t,CA)
fnplt(pp)
title('Concentration vs time')
legend('Experiment','Fitting')
dp=fnder(pp)
dCAdt=-fnval(dp,t);
y=log(dCAdt);x=log(CA);
p=polyfit(x,y,1);
beta0=[exp(p(2)),p(1)];
rate=@(beta,C)(beta(1)*C.^beta(2))
beta=nlinfit(CA,dCAdt,rate,beta0)
figure
plot(t,dCAdt,'o',t,rate(beta,CA),'*',t,rate(beta0,CA),'p')
title('Reaction rate vs time')
legend('Rate from difference','Predicted by nonlinear','Predicted by linear fitting')

从图中可以看出，三次样条较好拟合了原始实验数据，微分计算得到的反映速率存在拐点，非线性拟合在高反映速率时拟合效果好于线性拟合效果。

数值积分算法

积分是工程领域中常见数学计算。在微积分中,计算连续函数f(.x)在区间[a,b]上的积分是通过解析法求f(x)原函数F(x)得到的。

但在有些情况下,无法得到积分的解析解。如:

(1)被积函数以一组数据形式表示;

(2）被积函数过于特殊或原函数不能用初等函数表示,积分表中无法找到可沿用的现成公式;

(3)有的原函数十分复杂难以计算。

这时,用牛顿-莱布尼茨公式求定积分会失效,需要借助于数值积分法去解决问题。常用数值积分基本思路来自于插值法,它通过构造一个插值多项式P(x)作为f(x)的近似表达式,用P(x)的积分值作为f(x)近似积分值。数值积分的方法很丰富,常用插值型求积公式有两类:一类是等距节点的牛顿-柯特斯求积公式;另一类是不等距节点的高斯型求积公式。

一维数值积分函数包括: trapz和函数名以quad函数,integral(二维和三维数值积分函数为integral2和 integral3)。

1、复合梯形法数值积分函数: trapz调用形式:Z=trapz(X,Y)

输入变量X表示自变量,Y表示因变量的值;缺省参数X时,表示X被等分,步长为l;Z代表返回的积分值。

2．自适应辛普森法数值积分函数:quad基本调用格式:q=quad(fun，a，b,tol,trace,p1,p2)

fun——被积函数。可以是匿名函数,内联函数或函数句柄。无论采用哪种形式,函数应能接受一个向量x,并返回对应的函数值向量y。当采用函数句柄时,函数的声明语句应具有如下形式:

function y=fun(x，pl，p2，...)

输入变量x为自变量，pl、p2等为quad传递到fun中的参数;返回值y为积分表达式在x处的函数值。

a，b分别是积分的下限和上限。q—积分结果。

tol———-默认误差限﹐默认值为1.e—6。

trace——一取О表示不用图形显示积分过程,非0表示用图形显示积分过程。

pl，p2，...—-直接传递给函数fun的参数。

用trapz函数和quad函数求积分

 %使用trapz函数求积分

x=2:0.15:5;
y=log(x)./x.^2;
l1=trapz(x,y)
%使用quad函数求积分
f=@(x)(log(x)./x.^2)
I2=quad(f,2,5)
l1 =
    0.3247
f =
  包含以下值的 function_handle:
    @(x)(log(x)./x.^2)
I2 =
    0.3247

    可以看出使用trapz和quad函数求积分非常便捷，过冷水早期求积分的时候采用梯形法求积分，那个过程是相当麻烦，本期过冷水想要和大家分享的关于数值运算中微分和积分运算就这么多内容，希望上述知识程序能够有用，后期有遇到需要灵活处理的内容在做补充