MATLAB 数值微积分

---

前言

b站课程《MATLAB教程_台大郭彦甫（14课）》学习记录

---

一、Differentiation 微分

二、Polynomial Differentiation 多项式微分

三、Representing Polynomials in MATLAB

Polynomials were represented as row vectors 多项式用行向量表示
For example, consider the equation

To enter this polynomial into MATLAB, use
p = [1 0 -2 -5];

四、Values of Polynomials: polyval()

a = [9,-5,3,7]; x = -2:0.01:5;
f = polyval(a,x); 
plot(x,f,'LineWidth', 2);
xlabel('x'); ylabel('f(x)');
set(gca, 'FontSize', 14)

五、Polynomial Differentiation 多项式微分: polyder()

>> p=[5 0 -2 0 1];
polyder(p)
ans =
    20     0    -4     0
polyval(polyder(p),7)
ans =
        6832

六、Exercise

conv(): 多项式相乘

a = [20 -7 5 10];
b = [4 12 -3];
t = conv(a,b);
hold on;
x = -2:0.01:1;
plot(x,polyval(t,x),'--b');
plot(x,polyval(polyder(t),x),'-r');
hold off;
xlabel('x');
legend('f(x)','f`(x)');
set(gca,'XLim',[-2,1]);
set(gca,'YLim',[-800,800]);

七、Polynomial Integration 多项式积分

八、Polynomial Integration: polyint()

我们要给MATLAB提供一个常量3，作为C的值

>> 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

九、Numerical Differentiation 数值微分

十、Differences: diff()

diff()计算vector中相邻元素之间的差值

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

Exercise

obtain the slope（斜率） of a line between 2 points (1,5) and (2,7)

>> x = [1 2]; y = [5 7];
slope = diff(y)./diff(x)
slope =
     2

十一、Numerical Differentiation Using diff()

>> x0 = pi/2; h = 0.1;
x = [x0 x0+h];
y = [sin(x0) sin(x0+h)]; 
m = diff(y)./diff(x)
m =
   -0.0500

How does ℎ affect accuracy?
h越小越准确

Exercise

x0 = pi/2; h = 0.1;
while h>0.0000001
    x = [x0 x0+h];
    y = [sin(x0) sin(x0+h)]; 
    m = diff(y)./diff(x);
    h = h/10;
    m
end
m =
     0
m =
   -0.0500
m =
   -0.0050
m =
  -5.0000e-04
m =
  -5.0000e-05
m =
  -5.0000e-06
m =
  -5.0004e-07
m =
  -4.9960e-08

十二、How to Find the 𝑓′ over An Interval [0,2𝜋]?

十三、Find sin′(𝑥) over 𝑥 = [0, 2𝜋]

h = 0.5; x = 0:h:2*pi;
y = sin(x); m = diff(y)./diff(x);

十四、Various Step Size

The derivatives of 𝑓(𝑥) = sin(𝑥) calculated using various ℎ values

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, 'YLim', [0, 1.2]);
set(gca, 'FontSize', 18); 
%set(gca, 'FontName', 'symbol');
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;

Exercise

Given 𝑓(𝑥) = 𝑒 −𝑥 sin(𝑥 2/2), plot the approximate derivatives 𝑓′ of ℎ = 0.1, 0.01, and 0.001

g = colormap(lines); 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, 'Color', g(i,:));
end
hold off;
set(gca, 'XLim', [0, 2*pi]); set(gca, 'YLim', [-0.3, 0.3]);
set(gca, 'FontSize', 18);
set(gca, 'XTick', 0: 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;

十五、Second and Third Derivatives

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);

十六、Numerical Integration 数值积分

十七、Numerical Quadrature Rules 数值求积规则

十八、Midpoint Rule

十九、Midpoint Rule Using sum()

x(1:end-1) = [x1 x2 x3 x4 … xend-1]
x(2:end) = [x2 x3 x4 x5 … xend]

>> 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 =
   15.9950

二十、Trapezoid Rule

二十一、Trapezoid Rule Using trapz()

>> 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 =
   15.9950

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)
s =
   16.0100

二十二、Second-order Rule: 1/3 Simpson’s

二十三、Simpson’s Rule

>> 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 =
    16

二十四、Comparison

二十五、Review of Function Handles (@)

相当于嵌套函数
A handle is “a pointer to a function”
Can be used to pass functions to other functions
Example:
Pass a function f(x)=sin(x)
to a user-defined function: g(f,…)
f=sin(x)
g(@f,…)

二十六、Function Handles (@) Example

The input of the following function xy_plot is a math function:

function [y] = xy_plot(input,x)
% xy_plot receives the handle of a function
% and plots that function of x
y = input(x); plot(x,y,'r--');
xlabel('x'); ylabel('function(x)');
end

Try:

xy_plot(@sin,0:0.01:2*pi);

xy_plot(@cos,0:0.01:2*pi);

xy_plot(@exp,0:0.01:2*pi);

二十七、Numerical Integration: integral()

Numerical integration on a function from using global adaptive quadrature and default error tolerances
利用全局自适应正交和默认容错对函数进行数值积分

>> y = @(x) 1./(x.^3-2*x-5);
integral(y,0,2)
ans =
   -0.4605

二十八、Double and Triple Integrals

>> f = @(x,y) y.*sin(x)+x.*cos(y);
integral2(f,pi,2*pi,0,pi)
ans =
   -9.8696

>> f = @(x,y,z) y.*sin(x)+z.*cos(y);
integral3(f,0,pi,0,1,-1,1)
ans =
    2.0000

---

总结

继续加油吧。