微分

计算双曲正割函数的一阶导数和二阶导数

一阶导数和二阶导数的表达式是我们已知的，如下

首先定义计算区间和步长

1
2
dx = 0.1;   % 定义步长
x = -10:dx:10;  % 定义计算区间

按照解析式在 MATLAB 中绘制精确解

1
2
3
4
5
6
7
8
u = sech(x);
ux_exact = -sech(x).*tanh(x);
uxx_exact = sech(x)-2*sech(x).^3;

% 绘图
figure(1);
plot(x,u,x,ux_exact,x,uxx_exact);grid on
legend('u','ux','uxx');

双曲正割函数

要计算数值解，我们可以使用四阶的中心差分公式，在边界处使用向前差分公式和向后差分公式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
n  = length(x);
ux = zeros(1,n);

% 计算边界点
ux(1) = (-3*u(1)+4*u(2)-u(3))/(2*dx);
ux(2) = (-3*u(2)+4*u(3)-u(4))/(2*dx);

ux(n-1) = (3*u(n-1)-4*u(n-2)+u(n-3))/(2*dx);
ux(n) = (3*u(n)-4*u(n-1)+u(n-2))/(2*dx);

for j = 3:n-2
    ux(j) = (-u(j+2)+8*u(j+1)-8*u(j-1)+u(j-2))/(12*dx);
end

% 绘图
figure(2);
plot(x,ux_exact,x,ux,'o');grid on
legend('精确解','数值解');

双曲正割函数

计算二阶导数

1
2
3
4
5
6
7
8
9
% 绘制边界点
uxx(1) = (2*u(1)-5*u(2)+4*u(3)-u(4))/(dx*dx);
uxx(2) = (2*u(2)-5*u(3)+4*u(4)-u(5))/(dx*dx);
uxx(n-1) = (-2*u(n-1)+5*u(n-2)-4*u(n-3)+u(n-4))/(dx*dx);
uxx(n) = (-2*u(n)+5*u(n-1)-4*u(n-2)+u(n-3))/(dx*dx);

for j = 3:n-2
    uxx(j) = (-u(j+2)+16*u(j+1)-30*u(j)+16*u(j-1)-u(j-2))/(12*dx*dx);
end

数值微分

我们可以查看数值解的误差

一阶导数值解误差

二阶导数值解误差

误差精度均在 $10^{-5}$.

积分

积分程序我们可以直接使用 MATLAB 的内嵌函数，例如要使用梯形积分法则，就可以使用 trapz 函数，例如我们要求解 $u=\mathrm{sech}^2x$ 的积分，用一行代码就可以实现，结果为 2.0000 .

1
2
3
4
5
6
>> trapz(x,u.^2)

ans =

    2.0000

这个命令的计算误差在 $10^{-7}$ ，若要产生累积值，则可使用 cumtrapz 函数

1
2
3
4
int_sech = cumtrapz(x,u.^2);

figure(3);
plot(x,int_sech_sq);grid on

计算结果如图

sech(x)积分

如果我们需要计算一个指定区间内的积分，则可以使用 quad 函数，语法如下

1
int_quad = quad(@fun,a,b);

其中 @fun 是被积函数，积分下限为 a ，积分上限为 b . 直接调用是可行的，我们这里给出一个例子

1
2
A = 1;
int_quad = quad(inline('A*sech(x).^2','x','A'),-10,10,[],[],A);

这个例子没有直接调用，而是使用了 inline 函数，对参数进行了检查.

如果要计算二重积分或者三重积分，则可以使用 dblquad 函数.