学习笔记:郭彥甫 (Yan-Fu Kuo), 台大生機系 , MATLAB教學 - 10數值微積分
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()
- Plot the polynomial
- for .
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()
- Given
- What is the derivative(导数) of the function ?
- What is the derivative of the function value of ?
>> p = [5 0 -2 0 1];
>> polyder(p)
ans = 20 0 -4 0
>> polyval(polyder(p),7)
ans = 6832
Exercise
- Plot the polynomial and its derivative for .
Hint:
conv()
p1 = [20 -7 5 10]; p2 = [4 12 -3]; a = conv(p1,p2); x=-2:0.01:1; hold on plot(x,f,'b--','linewidth',2); plot(x,polyval(polyder(a),x),'r','linewidth',2); hold off % plot(x,f,'b--',x,polyval(polyder(a),x),'r','linewidth',2); xlabel('x') ylabel('y') legend('f(x)',"f'(x)");
Polynomial Intergration(积分):polyint()
for a polynomial
Given
What is the integral of the function with a constant of 3?
>> p = [5 0 -2 0 1]; >> polyint(p,3) ans = 1.0000 0 -0.6667 0 1.0000 3.0000
What is the derivative of the function value of ?
>> polyval(polyint(p,3),7)) ans = 1.6588e+04
Numerical Differentiation (数值微分):diff()
The simplest method: finite difference approximation(有限差分近似)
Calculating a secant line(割线) in the vicinty(附近) of
diff() : calculates the differences between adjacent(附近的) elements of a 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)
Given , find at using
>> x0 = pi/2 x0 = 1.5708 >> h=0.1 h = 0.1000 >> x = [x0 x0+h] x = 1.5708 1.6708 >> y = [sin(x0) sin(x0+h)] y = 1.0000 0.9950 >> m = diff(y)./diff(x) m = -0.0500
Exercise
Given , write a script to find the error of at using various
We have already known
x0 = pi/2; for h=[0.1 0.01 0.001 0.0001 0.00001 0.000001] x=[x0,x0+h]; y=[sin(x0),sin(x0+h)]; slop=diff(y)./diff(x); W=['h=',num2str(h),' error=',num2str(slop)]; disp(W) end
How to Find the over An Interval(区间) ?
- In the previous example, .
Strategy:
- Create an array in the interval
- The step is the .
- Calculate the at these points
For example:
h = 0.5; x = 0:h:2*pi;
Find over .
h = 0.5; x = 0:h:2*pi; y = sin(x); m = diff(y)./diff(x); % m = f'(x)
Various Step Size
The derivatives of caculated 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','p/4','p/2'}); % xtick是刻度(小竖线);xticklabel 刻度值(竖线下面的数值) h = legend('h=0.1','h=0.01','h=0.001','h=0.0001'); set(h,'FontName','Times New Roman'); box on;
Exercise
Given ,plot the approximate derivatives of , and
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,pi*2]); set(gca,'YLim',[-0.3,0.3]); set(gca,'FontSize',18); set(gca,'FontName','symbol'); set(gca,'XTick',0:pi/2:pi*2); set(gca,'XTickLabel',{'0','pi/2','pi','3/2pi','2pi'}); set(gca,'YTick',-0.2:0.1:0.2); set(gca,'YTickLabel',{'-0.2','-0.1','0','0.1','0.2'}); h = legend('h=0.1','h=0.01','h=0.001'); set(h,'FontName','Times New Roman'); box on;
Second and Third Derivatives(二阶、三阶导)
The second derivative and third derivative can be obtained using similar approaches
Given , plot and for .
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(数值积分)
Calculating the numerical value of a definite integral
Quadrature method — approximating the integral by using a finite set of points(正交方法 - 使用有限的点集近似积分)
Numerical Quadrature Rules(数值求积分准则)
Basic quadrature rules:
Midpoint rule(zeroth-order approximation)(中点准则(零阶近似),使用矩形预测面积大小)
Trapezoid rule(first-order approximation)(梯形准则(一阶近似),使用梯形预测面积大小)
Midpoint Rule Using sum()
Example:
>> h = 0.05; >> x = 0:h:2; % x is a vector >> midpoint = (x(1:end-1)+x(2:end))./2; % midpoint is a vector >> y = 4*midpoint.^3; % y is a vector >> s = sum(h*y) s =15.9950
Trapezoid Rule Using trapz()
Example:
>> h = 0.05; >> x = 0:h:2; >> y = 4*x.^3; >> s = h*trapz(y) 16.0100 % the forth step above equals trapezoid = (y(1:end-1)+y(2:end))/2; s = h*sum(trapezoid)
Second-order Rule:1/3 Simpson’s
Example:
>> 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
- Midpoint rule — zeroth-order
- Trapezoid rule — first-order
- Simpson’s rule — second-order
- higher order of approximation results to a higher accuracy (逼近的次数越高越精准)
Review of Function Handles(@)
- A handle is a pointer to a function
Can be used to pass function to other functions
For example, the input of the following function is another 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
When we need to use the function
xy_plot
,usexy_plot(@sin,0:0.01:2*pi);
Numerical Intergration:integral()
Numerical integration on a function from using global adaptive quadrature and default error tolerances.
Example:
>> y = @(x) 1./(x.^3-2*x-5); >> integral(y,0.2) ans = -0.4605
y = @(x) sin(x); integral(y,0,2)
Double and Triple Integrals(双重和三重积分)
Example:
f = @(x,y) y.*sin(x)+x.*cos(y); integral2(f,pi,2*pi,0,pi)
Example:
f = @(x,y,z) y.*sin(x)+z.*cos(y) integral3(f,0,pi,0,1,-1,1)