0
点赞
收藏
分享

微信扫一扫

MATLAB 数值微积分

西特张 2022-03-25 阅读 60

文章目录


前言

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

总结

继续加油吧。

举报

相关推荐

0 条评论