MATLAB台大郭彦甫老师课程笔记：第九课：数值微积分

1.MATLAB中表示多项式的方法：使用 Row vector

a = [9, -5, 3, 7];  x= -2: 0.01: 5 ;
f = polyval(a, x);    %% 表示计算在polynomial为a 的情况下，在x 处的数值
plot(x, f, 'LineWidth', 2);
xlabel('x'); ylabel('f(x)');
set(gca, 'FontSize', 14)

2.如何求一个多项式的微分：

>> p = [5 0 -2 0 1]; polyder(p)

ans =

    20     0    -4     0

>> polyval(polyder(p), 7)

ans =

        6832

第九课课后练习1：

%% 计算该polynomial 的值，及其微分数值
%% 提示：conv 是多项式相乘函数
a = [5 -7 5 10];  b = [4 12 -3];
c = conv(a, b);  x = -2: 0.01: 1;
f = polyval(c, x);
g = polyval(polyder(c), x);
hold on
plot(x, f, '--b', 'LineWidth', 2)
plot(x, g, '-.r', 'LineWidth', 2);
hold off
legend('f(x)', 'f''(x)'); %% 这里要打出f'(x)的话，只能用两个''，转置一下。
xlabel('x'); ylabel('f(x)');
set(gca, 'FontSize', 14)

那么，如何计算一个多项式的积分呢？
3.这时需要手动指定常数值（constant）。

>> p = [5 0 -2 0 1];
polyint(p, 3)

ans =

    1.0000         0   -0.6667         0    1.0000    3.0000

>> polyval(polyint(p, 3), 7)

ans =

   1.6588e+04

即：计算某个向量前后元素之间的差异。

>> x = [1 2 5 2 1];
diff(x)

ans =

     1     3    -3    -1

（这个diff 的函数可以用来算两点之间的斜率）

第九课课后练习2：
（更进一步）
How to find sin’(x) over x= [0, 2*pi]:

h = 0.5; x = 0: h: 2*pi;
y = sin(x); m = diff(y) ./ diff(x);
%% 由于这两条语句求出来的点数比原来的点数少最后一个点，因此可以将 m ：
% m = [diff(y) ./ diff(x) 1];  相当于手动添加最后一个点

diff指令是知道某点的函数值时候求导 polyder是知道多项式后求导

g = colormap(lines); hold on;
for i = 1:4
    x = 0: power(10, -i): pi;
    y = sin(x); m = diff(y)./ diff(x);
    plot(x(1: end-1), m, 'Color', g(i,:));
end
hold off;
set(gca, 'XLim', [0, pi/2]); 
set(gca, 'FontSize', 18);
set(gca, 'FontName', 'LaTex');
set(gca, 'XTick', 0:pi/4:pi/2);
set(gca, 'XTickLabel', {'0', '\pi/4', '\pi/2'});
h = legend('h=0.1', 'h=0.01', 'h=0.001', 'h=0.0001');
set(h, 'FontName', 'Times New Roman');
box on;

第九课课后练习3：

hold on;
for i = 1:3
    x = 0: power(10, -i): 2*pi;
    y = (exp(-x)).*sin((x.^2)/2); m = diff(y)./ diff(x);
    plot(x(1: end-1), m);
end
hold off;
set(gca, 'XLim', [0, 2*pi]); 
set(gca, 'FontSize', 18);
set(gca, 'FontName', 'LaTex');
set(gca, 'XTick', 0: pi/2: pi: 3*pi/2: 2*pi);
set(gca, 'XTickLabel', {'0', '\pi/2', '\pi', '3\pi/2', '2\pi'});
h = legend('h=0.1', 'h=0.01', 'h=0.001');
set(h, 'FontName', 'Times New Roman');
box on;

x = -2: 0.005: 2;  y = x .^3;
m = diff(y)./ diff(x);
m2 = diff(m)./diff(x(1: end-1));

plot(x, y, x(1: end-1), m, x(1: end-2), m2);
xlabel('x', 'FontSize', 18);
ylabel('y', 'FontSize', 18);
legend('f(x) = x^3', 'f''(x)', 'f''''(x)');
set(gca, 'FontSize', 18);

4.求积分的两种方法：

法一：中点法

实例用法：

h = 0.05; x = 0: h: 2;
midpoint = (x(1: end-1) + x(2: end)) ./2;
y = 4* midpoint .^3;
s = sum(h * y);

>> s

s =

   15.9950

Q: How to improve the accuracy?
A: 不断减小h 的值。

法二：梯形法
这个图的下面图上的公式好像有点问题，就只看上面的公式就好。

上题的实例用法：

h = 0.05; x = 0: h: 2;
y = 4 * x .^ 3;
s = h * trapz(y);

>> s

s =

   16.0100

%%Alternative:

h = 0.05; x = 0: h: 2;
y = 4 * x .^ 3;
trapezoid = (y(1: end-1) + y(2: end)) / 2;
s = h * sum(trapezoid)
% 梯形的面积公式：(上底+下底) * 高 / 2

5.相较之下更加精准的公式：辛普森公式

h = 0.05; x = 0: h: 2; y = 4*x.^3;
s = h/3*(y(1) + 2*sum(y(3: 2: end-2))+...
    4*sum(y(2: 2: end)) + y(end))

>> s
        
s =

    16

6.函数指针：

函数指针在计算式子的积分时的实例用法：

6.两变量和三变量积分：