【第3章】微積分問題的計算機求解

極限問題的解析解

單變量函數的極限

假設已知函數f(x),則極限問題一般描述為

\[L=\lim_{x \to x_0} f(x) \]

此外，還有單邊極限問題

\[L=\lim_{x \to x_0^-} f(x) \quad \text{左極限} \\ L=\lim_{x \to x_0^+} f(x) \quad \text{右極限} \]

matlab同樣可以做這些極限運算

L=limit(fun,x,x0)    %求極限
L=limit(fun,x,x0,'left'或'right')    %求單邊極限

舉個例子

試求解極限問題

\[\lim_{x \to \infty} \frac{sin(x)}{x} \]

matlab代碼

syms x; f=sin(x)/x; 
limit(f,x,0)

運行結果截圖

多變量函數的極限

若想求出二元函數的極限

\[L=\lim_{x \to x_0 \atop y \to y_0} f(x) \]

我們可以嵌套使用limit()函數。

L1=limit(limit(f,x,x0),y,y0) 或
L1=limit(limit(f,y,y0),x,x0)

如果x0或y0不是確定的值，而是另一個變量的函數，例如 $ x \to g(y) $ ，則上述極限求取順序不能交換。

函數導數的解析解

函數的導數和高階導數

如果函數fun和自變量x都已知，且均為符號變量，則可以用diff()函數解出給定函數的各階導數。

y=diff(fun,x)    %求導
y=diff(fun,x,n)  %求n階導

例：給出如下函數，試求出其一階導數

\[f(x)= \frac{sinx}{x^2+4x+3} \]

matlab代碼

syms x;
f=sin(x)/(x^2+4*x+3);
f1=diff(f,x,1);
latex(f1)

最后得出結果如下

\[\frac{\cos\!\left(x\right)}{x^2 + 4\, x + 3} - \frac{\sin\!\left(x\right)\, \left(2\, x + 4\right)}{{\left(x^2 + 4\, x + 3\right)}^2} \]

復合泛函求導

例：給出如下函數，試求出其三階導數公式，以及 $f(t)= e^{-t} $ 時的結果 關鍵 將f(t)聲明為符號變量

\[F(t) = t^2 sint f(t) \]

matlab代碼

syms t f(t); 
G=simplify(diff(t^2*sin(t)*f,t,3))
simplify(subs(G,f,exp(-t))), 
simplify(diff(t^2*sin(t)*exp(-t),3)-ans)

最后得出結果如下

G(t)=
6*cos(t)*f(t) + 6*sin(t)*diff(f(t), t) + 3*t^2*cos(t)*diff(f(t), t, t) + 12*t*cos(t)*diff(f(t), t) - 6*t*f(t)*sin(t) + t^2*sin(t)*diff(f(t), t, t, t) - 3*t^2*sin(t)*diff(f(t), t) + 6*t*sin(t)*diff(f(t), t, t) - t^2*cos(t)*f(t)

ans(t) =
2*exp(-t)*(3*cos(t) - 3*sin(t) + t^2*cos(t) + t^2*sin(t) - 6*t*cos(t))

ans(t) =
0

矩陣函數的求導

矩陣的求導：

\[H(x)= \begin{bmatrix} 4sin(5x) & e^{-4x^2} \\ 3x^2+4x+1 & \sqrt{4x^2+2} \\ \end{bmatrix} \]

可以對每個元素分別求導

syms x; 
H=[4*sin(5*x), exp(-4*x^2); 3*x^2+4*x+1, sqrt(4*x^2+2)];
H1=diff(H,x,3)

運行結果：

H1 =
 
[ -500*cos(5*x),                               192*x*exp(-4*x^2) - 512*x^3*exp(-4*x^2)]
[             0, (24*2^(1/2)*x^3)/(2*x^2 + 1)^(5/2) - (12*2^(1/2)*x)/(2*x^2 + 1)^(3/2)]

參數方程的導數

matlab中沒有直接能夠求解參數方程的函數，但我們可以根據其在數學上的定義來求：

根據遞推公式，我們可以從中看出了，它的形式與我們之前學習的遞歸調用有很大的相似性，因此我們可以編寫一個這樣的函數paradiff(y,x,t,n)來求參數方程的n階導數

%paradiff.m
function result=paradiff(y,x,t,n)
if mod(n,1)~=0, error('n should positive integer, please correct')
else, if n==1, result=diff(y,t)/diff(x,t);
   else, result=diff(paradiff(y,x,t,n-1),t)/diff(x,t);
end, end

例：已知參數方程如下，求其三階導數

\[\begin{cases} y=\frac{sint}{(t+1)^3} \\ x=\frac{cost}{(t+1)^3} \end{cases} \]

matlab代碼

syms t; 
y=sin(t)/(t+1)^3; 
x=cos(t)/(t+1)^3; 
f=paradiff(y,x,t,3); 
[n,d]=numden(f); %提取分子和分母，進行單獨化簡
F=simplify(n)/simplify(d)

運行結果

F =
(3*(t + 1)^7*(23*cos(t) + 24*sin(t) - 6*t^2*cos(t) - 4*t^3*cos(t) - t^4*cos(t) + 12*t^2*sin(t) + 4*t^3*sin(t) - 4*t*cos(t) + 32*t*sin(t)))/(3*cos(t) + sin(t) + t*sin(t))^5

多元函數的偏導數

matlab中沒有求取偏導數的專門函數，但我們仍可以通過diff()函數直接實現。假設已知二元函數f(x,y),若想求

\[\frac{\partial ^{m+n}}{\partial x^m \partial y^n}f(x,y) \]

則可以使用下面的函數求出

f=diff(diff(f,x,m),y,n) %或者
f=diff(diff(f,y,n),x,m)

例：求如下函數的兩個偏導數 $ \partial z / \partial x , \partial z / \partial y $

\[z=f(x,y)=(x^2-2x)e^{-x^2-y^2-xy} \]

matlab代碼

syms x y; 
z=(x^2-2*x)*exp(-x^2-y^2-x*y); 
zx=simplify(diff(z,x)), 
zy=diff(z,y)

運行結果

zx =
exp(- x^2 - x*y - y^2)*(2*x + 2*x*y - x^2*y + 4*x^2 - 2*x^3 - 2)
 
zy =
exp(- x^2 - x*y - y^2)*(- x^2 + 2*x)*(x + 2*y)

利用得到的偏導數函數zx，zy我們可以在z這個三維曲面上繪制出引力線，得到其梯度函數圖形表示，引力線物理意義可看作一個小球在這個位置所受的力。

[x0,y0]=meshgrid(-3:.2:2,-2:.2:2); 
z0=double(subs(z,{x,y},{x0,y0})); %將符號型轉為double型
surf(x0,y0,z0), zlim([-0.7 1.5]) %先畫出Z曲面
contour(x0,y0,z0,30), hold on; %畫出等高線並保持
zx0=subs(zx,{x,y},{x0,y0}); 
zy0=subs(zy,{x,y},{x0,y0}); %計算出各個點偏導數的值
quiver(x0,y0,-zx0,-zy0) %把偏導數結果用引力線形式表示出來

最終圖如下：

隱函數的偏導數

還是直接上結論吧，matlab沒有直接求隱函數的偏導數的函數，所以我們根據數學上的公式，編寫函數impldiff(f,x,y,n)對 $ z=f(x,y) $ 求n階偏導數

上代碼：

%impldiff.m
function dy=impldiff(f,x,y,n)
if mod(n,1)~=0, error('n should positive integer, please correct')
else, F1=-simplify(diff(f,x)/diff(f,y)); dy=F1;
   for i=2:n, dy=simplify(diff(dy,x)+diff(dy,y)*F1); 
end, end   

例：求如下二元隱函數的一階偏導數

\[z=f(x,y)=(x^2-2x)e^{-x^2-y^2-xy}=0 \]

matlab代碼

syms x y; 
f=(x^2-2*x)*exp(-x^2-y^2-x*y); 
F1=impldiff(f,x,y,1)

運行結果：

F1 =
(2*x + 2*x*y - x^2*y + 4*x^2 - 2*x^3 - 2)/(x*(x + 2*y)*(x - 2))

多元函數的雅可比(Jacobian)矩陣

假設有n個自變量的m個函數定義為

\[\begin{cases} y_1=f_1(x_1,x_2,\cdots x_n) \\ y_2=f_2(x_1,x_2,\cdots x_n) \\ \qquad \quad \vdots \qquad \vdots \\ y_m=f_m(x_1,x_2,\cdots x_n) \end{cases} \]

將相應的 $ y_i $ 對 $ x_i $ 求偏導，則得出矩陣

\[J= \begin{bmatrix} \frac{\partial y_1 }{\partial x_1 } & \frac{\partial y_1 }{\partial x_2 } & \cdots & \frac{\partial y_1 }{\partial x_n } \\ \frac{\partial y_2 }{\partial x_2 } & \frac{\partial y_2 }{\partial x_2 } & \cdots & \frac{\partial y_2 }{\partial x_n } \\ \vdots & \vdots & \ddots & \vdots \\ \frac{\partial y_m }{\partial x_1 } & \frac{\partial y_m }{\partial x_2 } & \cdots & \frac{\partial y_m }{\partial x_n } \end{bmatrix} \qquad \]

該矩陣又稱為雅可比(Jacobian)矩陣,在matlab中可以用jacobian()函數直接求得。該函數的調用格式為jacobian(x,y)，其中x是自變量構成的向量，y是由各個函數構成的向量。

例：已知直角坐標和極坐標變換公式如下：

\[x=rsin\theta cos\phi ,y=rsin\theta sin\phi, z=rcos\theta \]

求雅可比矩陣

matlab代碼

syms r theta phi; %三個自變量
x=r*sin(theta)*cos(phi); 
y=r*sin(theta)*sin(phi); 
z=r*cos(theta); %三個函數
J=jacobian([x; y; z],[r theta phi])

運行結果

J =
 
[ cos(phi)*sin(theta), r*cos(phi)*cos(theta), -r*sin(phi)*sin(theta)]
[ sin(phi)*sin(theta), r*cos(theta)*sin(phi),  r*cos(phi)*sin(theta)]
[          cos(theta),         -r*sin(theta),                      0]

積分運算

不定積分的求解

不定積分的形式

\[\int {f(x)}dx \quad \text{或者} \quad \int \cdots \int {f(x)}dx^n \]

matlab中用於計算積分的函數是int(fun,x)，其中fun為被積函數，x為積分變量。與求導數不同的是，當需要計算多重積分時，我們要使用嵌套的方式來計算多重積分。而且最終的到的結果是原函數，要在加常數C

\[F=int(\cdots int(fun,x)) \]

例：求如下函數的一階導數在對結果求積分，與原函數比較

\[f(x)= \frac{sinx}{x^2+4x+3} \]

matlab代碼

syms x; 
y=sin(x)/(x^2+4*x+3); 
y1=diff(y); 
y0=int(y1)

運行結果：

y0 =
 
sin(x)/(x^2 + 4*x + 3)

定積分與無窮積分的求解

\[\int_a^b {f(x)}dx \quad , \quad \int_a^\infty {f(x)dx} \]

matlab語句格式：

int(fun,x,a,b)
int(fun,x,a,inf)

還可以利用原函數計算

牛頓-萊布尼茲公式

\[\int_a^b {f(x)}dx =F(b)-F(a) \]

傅里葉(Fourier)級數逼近

給定周期函數 $ f(x), x \in[-L, L], T=2 L $ 可以人為的對該函數在其他區間上進行周期延拓，使得 $ f(x)=f(x+k T) $ ， k為任意整數，這樣可以根據需要將其寫成下面的級數形式

\[\begin{array}{c} f(x)=\frac{a_{0}}{2}+\sum_{n=1}^{\infty}\left(a_{n} \cos \frac{n \pi}{L} x+b_{n} \sin \frac{n \pi}{L} x\right) \end{array} \]

其中，

\[\left\{\begin{array}{ll} a_{n} & =\frac{1}{L} \int_{-L}^{L} f(x) \cos \frac{n \pi x}{L} \mathrm{~d} x, \quad n=0,1,2, \cdots \\ b_{n} & =\frac{1}{L} \int_{-L}^{L} f(x) \sin \frac{n \pi x}{L} \mathrm{~d} x, \quad n=1,2,3, \cdots \end{array}\right. \]

該級數稱為Fourier級數，而a_n，b_n又稱為Fourier系數。若x∈(a，b)，則可以計算出L=(b-a)/2，引入新變量金，使得x=x+L+a，則可以將f(x)映射成(-L，L)區間上的函數，可對之進行Fourier級數展開，再將金=x-L-a轉換成x的函數即可。MATLAB和Maple語言均未直提供求解Fourier系數與級數的現成函數。其實由上述公式不難編寫出解析或數值的Fourier級數求解函數。其中解析函數如下：

%求Fourier級數展開的MATLAB代碼 
%[F,A,B]=fseries(f,x,p,a,b)
function [F,A,B]=fseries(f,x,varargin)
[p,a,b]=default_vals({6,-pi,pi},varargin{:});
L=(b-a)/2; if a+b, f=subs(f,x,x+L+a); end
A=int(f,x,-L,L)/L; B=[]; F=A/2;
for n=1:p
   an=int(f*cos(n*pi*x/L),x,-L,L)/L; bn=int(f*sin(n*pi*x/L),x,-L,L)/L; 
    A=[A, an]; B=[B,bn]; F=F+an*cos(n*pi*x/L)+bn*sin(n*pi*x/L);
end
if a+b, F=subs(F,x,x-L-a); end

支持函數的編寫

function varargout=default_vals(vals,varargin)
if nargout~=length(vals), error('number of arguments mismatch'); 
else, nn=length(varargin)+1;
    varargout=varargin; for i=nn:nargout, varargout{i}=vals{i};
end, end, end

該函數的調用格式為〔A，B，F]=fseries(f，x，p，a，b)，其中，f為給定函數，x為自
變量，p為展開項數，a，b為x的區間，可以省略取其默認值[-π，π]，得出的A，B為
Fourier系數，F為展開式。該函數應該置於@sym路徑下。仿照解析解fseries(）函數
其實不難寫出其數值版，有興趣的讀者可以自己試一試。

例：試求給定函數：$ x(x-π)(x-2π),x\in (0,2π) $ 的傅里葉級數展開

matlab代碼：

syms x, y=x*(x-pi)*(x-2*pi);
F=fseries(y,x,12,0,2*pi);
F
%運行結果
F =

(3*sin(2*x))/2 + (4*sin(3*x))/9 + (3*sin(4*x))/16 + (12*sin(5*x))/125 + sin(6*x)/18 + (12*sin(7*x))/343 + (3*sin(8*x))/128 + (4*sin(9*x))/243 + (3*sin(10*x))/250 + (12*sin(11*x))/1331 + sin(12*x)/144 + 12*sin(x)

也就是這個：

\[f(x)= \sum _{n=1}^{\infty}\frac{12}{n^{3}}\sin nx \]

泰勒(Taylor)冪級數展開

單變量函數的泰勒冪級數展開

如果在x=0點處進行冪級數展開，那么有

\[f(x)=a_{1}+a_{2}x+a_{3}x^{2}+ \cdots +a_{k}x^{k-1}+o(x^{k}) \]

其中，系數$ a_i $可以由下面的公式求出來

\[a_i = \frac{1}{i!}\lim _{x \rightarrow 0}\frac{d^{i-1}}{dx^{i-1}}f(x),i=1,2,3\cdots \]

該冪級數展開又稱為Maclaurin級數，若關於x=a點進行展開，則可以得出

\[f((x)=b_{1}+b_{2}(x-a)+b_{3}(x-a)^{2}+ \cdots +b_{k}(x-a)^{k-1}+o \left[(x-a)^{k}\right] \]

\[b_{i}= \frac{1}{i!}\lim _{x \rightarrow a}\frac{d^{i-1}}{dx^{i-1}}f(x),i=1,2,3\cdots \]

matlab解法：

%在x=0點處展開
F=taylor(fun,x,'Order',k)
F=taylor(fun,x, k) %早期matlab版本
%在x=a點處展開
F=taylor(fun,x,a,'Order', k)
F=taylor(fun,x,k,a) %早期matlab版本

例：求下列函數在x=0,x=a處的展開

\[f(x)= \frac{sinx}{x^{2}+4x+3} \]

syms x; f=sin(x)/(x^2+4*x+3); 
y=taylor(f,x,'Order',9)
syms a; taylor(f,x,a,'Order',5)

結果太長了我就不寫了

多變量函數的泰勒冪級數展開

\[f(x)=f(a)+ \left[(x_{1}-a_{1})\frac{\partial}{\partial x_{1}}+ \cdots +(x_{n}-a_{n})\frac{\partial}{\partial x_{n}}\right] f(x)|_{x=a} \]

\[\frac{1}{2!}\left[(x_{1}-a_{1})\frac{\partial}{\partial x_{1}}+ \cdots +(x_{n}-a_{n})\frac{\partial}{\partial x_{n}}\right] ^{2}f(x)|_{x=a} + \cdots + \]

\[\frac{1}{k!}\left[(x_{1}-a_{1})\frac{\partial}{\partial x_{1}}+ \cdots +(x_{n}-a_{n})\frac{\partial}{\partial x_{n}}\right] ^{k}f(x)|_{x=a} + \cdots \]

matlab求法：

F=taylor(f,[x1,x2,···,xn],[a1,a2,···,an」,'Order',k)

例：二元函數

\[z=f(x,y)=(x^{2}-2x)e^{-x^{2}-y^{2}-xy} \]

1.在原點展開Taylor級數

symsxy；f=(x^2-2*x)*exp(-x^2-y^2-x*y)；
F=taylor(f,[x,y],[0,0],'Order’,8)；
collect(F,x)

2.關於(1,a)點展開

syms_a；
F=taylor(f,[x,y],[1,a],'Order’,3),
F1=simplify(F)

序列求和與序列求積運算