MATLAB（九）数值微积分

前言

此篇文章是我在B站学习时所做的笔记，主要讲多项式微分与积分和数值微分与积分，部分为亲自动手演示过的，方便复习用。此篇文章仅供学习参考。

---

提示：以下是本篇文章正文内容，下面案例可供参考

用MATLAB表示多项式

---

---

多项式的值：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)

---

---

多项式微分：polyder()

1. 函数f’(x)的导数是什么?

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

ans =

    20     0    -4     0

意思就是f’(x)=20*x³-0*x²-4*x¹+0*x⁰

2. 函数f’(7)的导数值是什么?

>> polyval(polyder(p),7)

ans =

        6832

---

---

Exercise练习

a=conv([5 -7 5 10],[0 4 12 -3]); 
%polyder(a); 
hold on
d=linspace(-2,1);
c=polyval(a,d);
b=polyval(polyder(a),d) 
e=plot(d,c);
f=plot(d,b); 
set(e,'LineWidth',2,'linestyle','--','color','b')
set(f,'color','r') 
g='f(x)';
h="f'(x)" ;
legend(g,h)
%set(gca, 'linewidth', 2)
box on

---

---

多项式的集成：polyint()

1. 对常数为3的函数∫f(x)积分是什么？

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

ans =

    1.0000         0   -0.6667         0    1.0000    3.0000
>> format rational
>> polyint(p, 3)
ans =
       1              0             -2/3            0              1              3  

意思是他的积分等于 1x⁵+0x⁴-2/3x³+0x²+1*x¹+3

2. ∫f(7)的函数值的导数是什么?

>> polyval(polyint(p, 3),7)
ans =
   49765/3    

---

---

偏差：diff()

• diff() 计算向量中相邻元素之间的差值

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

---

---

Exercise练习

求两点(1,5)和(2,7)之间直线的斜率

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

---

---

数值微分法使用 diff()

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

---

---

如何在区间[0,2m]中求f’ ?

---

---

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

x(1:end-1)就是去掉x中的最后一个值，保证画图时两个变量数目一致

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

---

---

不同的步长

f(x) = sin(x)的导数用不同的h值计算
power(10, -i)表示为10^(-i)

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, 'XTick', 0:pi/4:pi/2);
set(gca, 'XTickLabel', {'0', 'p/4', 'p/2'});
h = legend('h=0.1','h=0.01','h=0.001','h=0.0001');
set(h,'FontName', 'Times New Roman'); box on;

---

---

Exercise练习

---

---

二阶和三阶导数

• 二阶导数f’‘和三阶导数f’’'可以用类似的方法得到
• 已知f(x)=x³,作图使f’ 和f’'在-2≤x≤2

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)','Location',"southeast");
set(gca, 'FontSize', 18);

---

---

数值积分

• 计算定积分的数值
  
• 求积法-利用有限点集逼近积分
  

---

---

数值求积规则

基本正交规则：

1. 中点规则(零阶近似)
   
2. 梯形法则(一阶近似)
   

---

---

中点规则

---

---

使用sum()的中点规则

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

s = 15.9950

想让它计算更精准，就调整h的值，使它变得更小，如h=0.0001

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

s = 16.0000

---

---

梯形法则

---

---

使用trapz()的梯形规则

trapz() 梯形数值积分

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

s = 16.0100

或者另外一种方式

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

---

---

二阶规则：1/3Simpson’s

---

---

辛普森法则

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

---

---

比较

当逼近的次数越高，越精确，误差越小。

---

---

函数句柄(@）

• 下面函数xy plot的输入是一个数学函数:

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

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

结果：

---

---

数值积分: integral()

利用全局自适应求积和默认误差公差对函数进行数值积分

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

---

---

二重积分和三重积分

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

---

如若侵权，请及时与我联系。