MATLAB教学_10数值微积分

最新推荐文章于 2022-07-13 21:55:38 发布

这么大个圈圈

最新推荐文章于 2022-07-13 21:55:38 发布

阅读量2.7k

点赞数 13

分类专栏： MATLAB

本文链接：https://blog.csdn.net/bsjsosososo/article/details/104471837

版权

MATLAB 专栏收录该内容

12 篇文章 95 订阅

订阅专栏

本文学习视频地址：https://www.bilibili.com/video/av68228488?p=10

课堂PPT以及本人学习代码已上传。

本文学习内容：

多项式的微分和积分
数值的微分和积分

Numerical Differentiation

多项式的表示方法

例如： f(x)=x^3-2x-5;

也就是 x^3 的系数为1， x^2的系数为0， x 的系数为 -2，常数为-5；

所以 p=[1 0 -2 -5];

polyval()

y = polyval(p,x) 计算多项式 p 在 x 的每个点处的值。参数 p 是长度为 n+1 的向量，其元素是 n 次多项式的系数（降幂排序）：

p(x)=p1^xn+p2^(xn−1)+...+pn^x+pn+1.

虽然可以为不同目的使用 polyint、polyder 和 polyfit 等函数计算 p 中的多项式系数，但您也可以为系数指定任何向量。

要以矩阵方式计算多项式，请改用 polyvalm。

[y,delta] = polyval(p,x,S) 使用 polyfit 生成的可选输出结构体 S 来生成误差估计值。delta 是使用 p(x) 预测 x 处的未来观测值时的标准误差估计值。

clear all;   %6分钟练习
clc;
a=[9,-5,3,7];     %构建多项式的系数
x=-2:0.01:5;
f=polyval(a,x);   %计算多项式 f在每个x 点的值
plot(x,f,'LineWidth',2,'Color','red'); %用图表示出来
xlabel('x'); ylabel('f(x)');
set(gca,'FontSize',14);

polyder()

多项式微分

clear all;  %11分钟练习
clc;
p=[5 0 -2 0 1];
polyder(p)        %f'(p)
polyval(polyder(p),7)   % f'(7)

16分钟练习

做出来和老师的不一样。有的弹幕说是 p=[5 -7 5 10]; 这里想不通。

conv()

w = conv(u,v) 返回向量 u 和 v 的卷积。如果 u 和 v 是多项式系数的向量，对其卷积与将这两个多项式相乘等效。

w = conv(u,v,shape) 返回如 shape 指定的卷积的分段。例如，conv(u,v,'same') 仅返回与 u 等大小的卷积的中心部分，而 conv(u,v,'valid') 仅返回计算的没有补零边缘的卷积部分。

clear all;  %16分钟练习
clc;
p=[20 -7 5 10];
q=[4 12 -3];
x=-2:0.01:1;
f=conv(p,q);    %如果是两个相乘的是行向量，则可以直接使用 conv()
a=polyval(f,x);   %计算 f(x)的值
b=polyval(polyder(f),x);  %计算 f'(x)的值
plot(x,a,'--b',x,b,'red','LineWidth',2);
xlabel('x'); 
legend('f(x)',"(f'(x))")

polyint()

多项式积分

q = polyint(p,k) 使用积分常量 k 返回 p 中系数所表示的多项式积分。

q = polyint(p) 假定积分常量 k = 0。

clear all;      %22分钟练习
p=[5 0 -2 0 1];
x=-2:0.01:8
polyint(p,3)            %求积分，常数指定为3
f=polyval(polyint(p,3),x);   %计算积分的值
plot(x,f);      
polyval(polyint(p,3),7)    %求x=7时积分的值

Numerical Differentiation

diff()

差分和近似导数

Y = diff(X) 计算沿大小不等于 1 的第一个数组维度的 X 相邻元素之间的差分：

如果 X 是长度为 m 的向量，则 Y = diff(X) 返回长度为 m-1 的向量。Y 的元素是 X 相邻元素之间的差分。
```
Y = [X(2)-X(1) X(3)-X(2) ... X(m)-X(m-1)]
```
如果 X 是不为空的非向量 p×m 矩阵，则 Y = diff(X) 返回大小为 (p-1)×m 的矩阵，其元素是 X 的行之间的差分。
```
Y = [X(2,:)-X(1,:); X(3,:)-X(2,:); ... X(p,:)-X(p-1,:)]
```
如果 X 是 0×0 的空矩阵，则 Y = diff(X) 返回 0×0 的空矩阵。

Y = diff(X,n) 通过递归应用 diff(X) 运算符 n 次来计算第 n 个差分。在实际操作中，这表示 diff(X,2) 与 diff(diff(X)) 相同。

Y = diff(X,n,dim) 是沿 dim 指定的维计算的第 n 个差分。dim 输入是一个正整数标量。

clear all;      %30分钟练习
f=[1 2 5 2 1];   
diff(f)          %diff()返回的是 [2-1 5-2 2-5 1-2]
x=[1 2];        %因为要求斜率，所以根据坐标，将X和Y的坐标重新规划
y=[5 7];     
slope=diff(y)./diff(x)    %这里用点乘，以防x 和 y是多个数据的矩阵

clear all;      %34分钟练习
% x=0:pi/50:2*pi;
% f=sin(x);
% g=polyder(f);       %这里不能用polyder(),因为这个函数是多项式求导，并不是连续函数
% polyval(g,x)
x0=pi/2; h=0.1;
x=[x0 x0+h];
y=[sin(x0) sin(x0+h)];
m=diff(y)./diff(x)

39分钟练习

x0=pi/2; h=0.1;    %39分钟练习
for i=1:6
    x=[x0 x0+h];
    y=[sin(x0) sin(x0+h)];
    m(i)=diff(y)./diff(x);
    h=h/10;
end
format long
m

45分钟练习

%%
h=0.5;  x=0:h:2*pi;   %45分钟练习
y=sin(x);   
m=diff(y)./diff(x);
plot(x,y);
x(length(x))=[];   %弹幕大神，把x的最后一项设置为空。以此和plot()匹配
hold on
plot(x,m)

%%
h=0.5;  x=0:h:2*pi;    %另一个方法
y=sin(x);   
m=diff(y)./diff(x);
plot(x,y);
hold on
plot(x(:,1:end-1),m)   %这里用到end的一个用法：索引
%x(:,1:end-1), 冒号是表示所有行，因为在第一个属性里，第二个1：end-1 表示 从第一列到倒数第二列。

使用 end 访问矩阵 A 的最后一行。

A = magic(3)

A = 3×3

     8     1     6
     3     5     7
     4     9     2

B = A(end,1:end)      %这里第一个end 是最后一行， 第二个 1：end 是指列从第1列到最后一列

B = 1×3

     4     9     2

具体可以查看：https://ww2.mathworks.cn/help/matlab/matlab_oop/object-end-indexing.html

不过我没看懂。

clear all;    %46分钟练习
g=colormap(lines);    %首先是将内置的colormap(lines)赋值给了g, lines是里面内置的一种线的颜色矩阵
hold on;
for i=1:4
    x=0:power(10,-i):pi;    % power()是按元素尔幂函数，这里是10^(-i)
    y=sin(x);   m=diff(y)./diff(x);   
    plot(x(:,1:end-1),m,'Color',g(i,:));   %一共有4条线，依次取了lines里面第1行到第4行的颜色
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','\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;

C = power(A,B) 是执行 A.^B 的替代方法，但很少使用。它可以启用类的运算符重载。

47分钟练习

这里exp(x) ; 就是e^x;

clc;
clear all; %47分钟练习
g=colormap(flag);
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
plot(x,y,'Color',g(4,:));
hold off;
xlim([0 2*pi]);
set(gca,'FontSize',20);
set(gca,'XTick',0:pi/2:2*pi);
set(gca,'XTickLabel',{'0','\pi/2','\pi','3\pi/2','2\pi'});
legend('h=0.1','h=0.01','h=0.01','h=0.001','y(x)');
xlabel('x');
ylabel('y');
box on;

多次积分

思路是先算出一次的积分，再用一次的积分再求积分。

clear all; clc; %49分钟练习
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);

数值积分

区间积分也就是求面积。有两种方法：

Midpoint rule ：

就是取中间点。在这个方法里面 h，也就是宽度是一定的。只有高度不一样。高度也就是f(x)的值。

下面的图，老师讲了怎么获取这个中间点的办法。第1个 x 矩阵是从x1开始，第二个矩阵是从x2开始，这两个相加除2就是中间点。

clear all; clc; %59分钟练习
h=0.05; x=0:h:2;
midpoint=(x(1:end-1)+x(2:end))./2;   %求积分也就是求面积，算出每个中间点
y=4*midpoint.^3;   % 用中间点带入到积分中，得出y值；
s=sum(h*y)   %求出用宽度 h和高度 y的乘积并相加
        %如果要精度更高，可以把 h 取小

Trapezoid rule

是用梯形的办法求出来。

Q = trapz(Y) 通过梯形法计算 Y 的近似积分（采用单位间距）。Y 的大小确定求积分所沿用的维度：

如果 Y 为向量，则 trapz(Y) 是 Y 的近似积分。
如果 Y 为矩阵，则 trapz(Y) 对每列求积分并返回积分值的行向量。
如果 Y 为多维数组，则 trapz(Y) 对其大小不等于 1 的第一个维度求积分。该维度的大小变为 1，而其他维度的大小保持不变。

clear all; clc;   %62分钟练习
h=0.05; x=0:h:2;    
y=4*x.^3;
s=h*trapz(y)      %trapz()来计算梯形数值积分
%%
clear all; clc;   %64分钟练习
h=0.05; x=0:h:2;    
y=4*x.^3;
trapezoid=(y(1:end-1)+y(2:end))/2;
s=h*sum(trapezoid)

辛普森的方法会准确一点，公式不会。大家自行百度。

三种方法的比较：

Function Handles(@)

句柄就像是一个指向一个函数的指针，在使用的时候利用这个指针指向其它函数，就可以调用。

例如： f=sin(x); 如果我要在函数g() 里用 f时，那么就要用指针去指向sin(x). 如下： g(@f,...)

创建函数句柄

通过在函数名称前添加一个 @ 符号来为函数创建句柄。例如，如果您有一个名为 myfunction 的函数，请按如下所示创建一个名为 f 的句柄：

f = @myfunction;

使用句柄调用函数的方式与直接调用函数一样。例如，假设您有一个名为 computeSquare 的函数，该函数定义为：

function y = computeSquare(x)
y = x.^2;
end

创建句柄并调用该函数以计算 4 的平方。

f = @computeSquare;     %%相当于调用
a = 4;
b = f(a)

b =

    16

匿名函数

您可以创建指向匿名函数的句柄。匿名函数是基于单行表达式的 MATLAB 函数，不需要程序文件。构造指向匿名函数的句柄，方法是定义 anonymous_function 函数主体，以及指向匿名函数 arglist 的以逗号分隔的输入参数列表。语法为：

h = @(arglist) anonymous_function

也就是在（arglist）里写未知写，后面紧跟着写函数

例如，创建一个指向用于计算平方数的匿名函数的句柄 sqr，并使用其句柄调用该匿名函数。

sqr = @(n) n.^2;
x = sqr(3)

x =

     9

clear all; clc;   %75分钟练习
y=@(x) 1./(x.^3-2*x-5);     %%有一个输入数
integral(y,0,2)
f=@(x,y) y.*sin(x)+x.*cos(y);   %%有两个输入数
integral2(f,pi,2*pi,0,pi)
f=@(x,y,z) y.*sin(x)+z.*cos(y);  %%有三个输入数
integral3(f,0,pi,0,1,-1,1)

这么大个圈圈

关注

13
点赞
踩
41

收藏

觉得还不错? 一键收藏
2
评论
MATLAB教学_10数值微积分

本文学习视频地址：https://www.bilibili.com/video/av68228488?p=10课堂PPT以及本人学习代码已上传。本文学习内容：多项式的微分和积分数值的微分和积分目录多项式的表示方法polyval()polyder()16分钟练习conv()polyint()Numerical Differentiation39分钟练...
复制链接

扫一扫