學習matlab之微積分
🍎微積分部分
極限
理論
limit(f,x,a)=limit(f,a) %求x->a時的極限
limit(f) = limit(f,x,0) %默認趨於零
limit(f,x,a,'right') %右極限
limit(f,x,a,'left') %左極限
limit(f,x,inf) %無窮
實驗
>> syms x;
>> f=x*sin(1/x);
>> limit(f,x,0)
ans =
0
>> g=sin(x)/x;
>> limit(g)
ans =
1
>> x=-1:0.01:1;
>> y=x.*sin(1./x);
>> plot(x,y)
>> syms x; %下面這個實驗證明了左極限等於右極限，極限存在
>> y=x*sin(1/x);
>> limit(y,x,0,'right')
ans =
0
>> limit(y,x,0,'left')
ans =
0
>> limit(y)
ans =
0
>> syms x;
>> f=sin(1/x);
>> plot(x,y) %這個錯誤證明了不能顯示非數組類型的圖像，因此之前定義syms為系統變量
錯誤使用 plot
數據必須為可轉換為雙精度值的數值、日期時間、持續時間或數組。
>> limit(f)
ans =
NaN %不存在
模型舉例：
預測某地區耐用消費品的市場保有量：
假設： (1) 第 n 個月，商品的保有總量為 x(n) ；
(2) 每個月每個用戶會影響 λ 個人購買商品；
(3) λ 會隨着總數 x(n) 增加而減少；
(4) 不考慮其他因素影響.
我們不妨假設最簡單的關系，λ與x(n)是線性關系，可能之后我們要用數據集來修正模型，但即使這樣，我們也應該從簡單的入手：
\[x(n+1)=(1+\lambda)*x(n)\\ \lambda=a-b*x(n)\\ x(n+1)=(1+a-b*x(n))*x(n) \]
令y(n)=b*x(n),則有
\[y(n+1)=(1+a-y(n))*y(n) \]
老師說這是離散形式的阻滯增長模型,以此我們用matlab畫圖a的改變對圖像的影響，
function y=logistic(y0,a,n)
y=zeros(n,1);
y(1)=y0;
for i=2:1:n
y(i)=y(i-1)*(1+a-y(i-1));
end
plot(y)
導數
diff:直接求導
（我們將結合具體的例子來介紹基本求導操作）
\[求y=ln(x+cosx)的導數 \]
>> syms x;
>> diff(log(x+cos(x)))
ans =
-(sin(x) - 1)/(x + cos(x))
\[眾所周知，在matlab中，ln(x)寫為log(x) \]
diff:向量函數求導
在matlab中我們可以把幾個函數以向量的方式寫出，即定義一個向量函數，自然而然地，我們也可以對向量函數求導：
>> diff([sin(x)+1,cos(x)+x,4*x^6])
ans =
[ cos(x), 1 - sin(x), 24*x^5]
diff:高階導數
我們甚至還可以求高階導數，只需加個參數n即可
>> diff(log(x),4)
ans =
-6/x^4
diff:參數方程求導
也可以求帶參數方程形式的求導
\[x=x(t)\\ y=y(t)\\ \frac{dy}{dx}=\frac{y'(t)}{x'(t)} \]
\[x=t^2-ln(2+sin(t))\\ y=t^3-3*sin(ln(t))\\ 求\frac{dy}{dx}. \]
>> syms t;
>> dx_dt=diff(t^2-log(2+sin(t)));
>> dy_dt=diff(t^3-3*sin(log(t)));
>> dy_dx=dy_dt/dx_dt
dy_dx =
-((3*cos(log(t)))/t - 3*t^2)/(2*t - cos(t)/(sin(t) + 2))
極值和最值
fminbnd:區間上最值
找區間上的最小值,如果區間上有多個的話只返回一個點，具體返回的算法還沒探究過
>> [x,f]=fminbnd('sin(x)',0,10)
x =
4.7124
f =
-1.0000
>> [x,f]=fminbnd('sin(x)',-4,7)
x =
-1.5708
f =
-1.0000
>> [x,f]=fminbnd('sin(x)',0,30)
x =
10.9956
f =
-1.0000
fminsearch:某點附近最值
求函數在某一點附近的最小值也可以
>> [x,f]=fminsearch('sin(x)',pi/2)
x =
4.7124
f =
-1
但我們知道了matlab好像沒有專門找最大值的函數，因為我們只需給所求函數加個負號的工作沒有重新定義一個函數；
偏導數
diff:拿來驗證作業答案
一定要記住偏導數是針對多元函數而言的；
讓我們拿工科數學分析下冊里面的習題來練練手吧！
\[求u=\frac{1}{\sqrt{x^2+y^2+z^2}}的偏導數， \]
>> syms x;
>> syms y;
>> syms z;
>> du_dx=diff((x^2+y^2+z^2)^(-1/2),x)
du_dx =
-x/(x^2 + y^2 + z^2)^(3/2)
diff:高階偏導數
>> syms x y z;
>> diff(arctan(y/x),x,2)
函數或變量 'arctan' 無法識別。
是不是想輸入:
>> diff(atan(y/x),x,2)
ans =
(2*y)/(x^3*(y^2/x^2 + 1)) - (2*y^3)/(x^5*(y^2/x^2 + 1)^2)
經化簡上述答案可寫為(所以拿來驗證答案可能稍有不妥🥀)
\[\frac{2xy}{(x^2+y^2)^2} \]
diff:混合偏導數
>> diff(diff(atan(y/x),x),y)
ans =
(2*y^2)/(x^4*(y^2/x^2 + 1)^2) - 1/(x^2*(y^2/x^2 + 1))
>> diff(diff(atan(y/x),y),x)
ans =
(2*y^2)/(x^4*(y^2/x^2 + 1)^2) - 1/(x^2*(y^2/x^2 + 1))
再次證明了初等多元函數混合偏導數與順序無關，但其化簡做得確實不好（對於機器來說化簡是很難的吧，他們只會做有規律的事）
\[\frac{y^2-x^2}{(x^2+y^2)^2} \]
diff:多元向量函數偏導數
\[求u=\left ( \begin{array}{l} x^2+sin(y) \\ y^2+sin(z) \\ z^2+sin(x) \end{array} \right ) 的Jacobian矩陣 \]
>>syms x y z;
>> jacobian([x^2+sin(y),y^2+sin(z),z^2+sin(x)],[x,y,z])
ans =
[ 2*x, cos(y), 0]
[ 0, 2*y, cos(z)]
[ cos(x), 0, 2*z]
diff:隱函數形式函數求導
\[對於F(x,y)=0，我們可知\\ dF(x,y)=F_xdx+F_ydy=0\\ 即\frac{dy}{dx}=-\frac{F_x}{F_y}\\ 對於F(x,y,z)=0\\ \frac{\partial z }{\partial x} = -\frac{F_x}{F_z}\\ \]
>> clear
>> syms x y z;
>> F=x^2*exp(-2*y-2*z)-5;
>> dz_dy=-diff(F,y)/diff(F,z)
dz_dy =
-1
模型舉例
例5. 商場每 T 天購進 Q 公斤某種蔬菜零售，每天銷量為 r
公斤，rT > Q. 那么 T 和 Q 如何選擇才能讓總支出最小化？
假設：
(1) 每次進貨需要額外開支a 元；
(2) 每天每公斤蔬菜儲存費支出b 元；
(3) 每天每公斤缺貨損失費c 元
\[我們不妨記q(t)為t時刻的蔬菜存量，具體形式：\\ q(t)=Q-rt,t\in[0,T],\\ 可知t=\frac{Q}{r}時蔬菜售完，則在之后的(\frac{Q}{r},T]缺貨，此時q(t)<0,\\ 則一個訂貨周期的總支出為：\\ \begin{equation*} \begin{aligned} C &=a+b*\int_{0}^{\frac{Q}{r}} q(t)dt+c*\int_{\frac{Q}{r}}^{T} |q(t)|dt\\ &=a+b*(Qt-\frac{1}{2}rt^2)|^{\frac{Q}{r}}_{0}-c*(Qt-\frac{1}{2}rt^2)|^{T}_{\frac{Q}{r}}\\ &=a+(b+c)*\frac{Q^2}{2r}-c*(QT-\frac{1}{2}rT^2)\\ &=a+b*\frac{Q^2}{2r}+\frac{c*(rT-Q)^2}{2r} \end{aligned} \end{equation*} \]
\[則每天的平均支出為C(T,Q)=\frac{a}{T}+\frac{b*Q^2}{2rT}+\frac{c*(rT-Q)^2}{2rT} \]
要使得C最小，則要使
\[\frac{\partial C(T,Q) }{\partial T}=0\\ \frac{\partial C(T,Q) }{\partial Q}=0 \]
>> syms a b c r T Q;
>> C=a/T+b*Q^2/(2*r*T)+c*(rT-Q)^2/(2*r*T);
函數或變量 'rT' 無法識別。
>> C=a/T+b*Q^2/(2*r*T)+c*(r*T-Q)^2/(2*r*T);
>> dC_dT=diff(C,T)
dC_dT =
- a/T^2 - (c*(Q - T*r))/T - (c*(Q - T*r)^2)/(2*T^2*r) - (Q^2*b)/(2*T^2*r)
>> dC_dQ=diff(C,Q)
dC_dQ =
(c*(2*Q - 2*T*r))/(2*T*r) + (Q*b)/(T*r)
>> [Q,T]=solve(dC_dT,dC_dQ,Q,T)
Q =
(2^(1/2)*c*r*((a*(b + c))/(b*c*r))^(1/2))/(b + c)
-(2^(1/2)*c*r*((a*(b + c))/(b*c*r))^(1/2))/(b + c)
T =
2^(1/2)*((a*(b + c))/(b*c*r))^(1/2)
-2^(1/2)*((a*(b + c))/(b*c*r))^(1/2)
手算解得
\[Q=\sqrt{\frac{2a}{b}\frac{rc}{b+c}}\\ T=\sqrt{\frac{2a}{b}\frac{b+c}{rc}} \]
按照這個來進貨可以使得每天平均支出最小化，即總支出最小化
導數的應用
驗證洛必達法則
\[求\lim_{x \to 0} \frac{3^x-2^x}{x} \]
>> syms x;
>> f=(3^x-2^x)/x;
>> limit(f)
ans =
log(3) - log(2)
>> f=(3^x-2^x);
>> g=x;
>> limit(diff(f,x)/diff(g,x),x,0)
ans =
log(3) - log(2)
驗證Taylor展開
taylor(f,x,'ExpansionPoint',a,'Order',n)
tricks:當matlab遇到函數忘記函數參數時可按CRTL+F1
>> clear
>> syms x f
>> f=exp(x);
>> T2=taylor(f,x,3)
T2 =
exp(3) + exp(3)*(x - 3) + (exp(3)*(x - 3)^2)/2 + (exp(3)*(x - 3)^3)/6 + (exp(3)*(x - 3)^4)/24 + (exp(3)*(x - 3)^5)/120
>> T2=taylor(f,x,0)
T2 =
x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1
>> clear
>> syms x;
>> f=exp(x);
>> T1=taylor(f,x,'ExpansionPoint',0,'Order',2);
>> T2=taylor(f,x,'ExpansionPoint',0,'Order',3);
>> T3=taylor(f,x,'ExpansionPoint',0,'Order',4);
>> plot(x,f,x,T1,x,T2,x,T3)
錯誤使用 plot
數據必須為可轉換為雙精度值的數值、日期時間、持續時間或數組。
>> plot(x,f,'^',x,T1,'*',x,T2,'-',x,T3)
錯誤使用 plot
數據必須為可轉換為雙精度值的數值、日期時間、持續時間或數組。 % 這里的plot只接受上面類型的函數，並不接受syms變量
>> fplot(f,'-') % 如果要用syms畫圖，就得用fplot,這樣的好處是精度夠高
>> fplot(T1)
>> fplot(f,'r')
>> hold on
>> fplot(T1,'b')
>> fplot(T2,'g')
>> fplot(T3)
>> t='$y=exp(x)$ and some of its taylor style ';
>> title(t,'interpreter','latex') % 顯示latex題目的方法
解方程
solve?
\[f(x)=0 \]
>> clear
>> syms x;
>> f=x^2+3*x+2;
>> solve(f) % 這樣是默認解f=0
ans =
-2
-1
\[\left\{\begin{matrix} f(x,y)=0\\ g(x,y)=0 \end{matrix}\right. \]
>> clear
>> syms x y a1 b1 c1 a2 b2 c2;
>> f=a1*x+b1*y+c1;
>> g=a2*x+b2*y+c2;
>> [x,y]=solve(f,g,x,y) %無論是符號解還是精確解都可以解出
x =
(b1*c2 - b2*c1)/(a1*b2 - a2*b1)
y =
-(a1*c2 - a2*c1)/(a1*b2 - a2*b1)
這里其實可以發現一個問題，solve不一定會輸出問題的所有解
>> syms x;
>> f=sin(x);
>> solve(f)
ans =
0
>> solve(f-1/2)
ans =
pi/6
(5*pi)/6
>> solve(f-1)
ans =
pi/2
牛頓迭代法
牛頓迭代法是一個很重要的求解方程近似解的方法：
\[如果f(x)在[a,b]上二階可導，f(a)f(b)<0且f'(x) 與f''(x)在 在[a,b]上不變號，則可用牛頓迭代法來求解f(x)=0. \]
\[\begin{align*} &牛頓迭代法的具體步驟：\\ &(1) 輸入精度指標 ε>0;\\ &(2) 確定區間[a,b] ，滿足f(a)f(b)<0 且f'(x) 與f''(x) 不變號;\\ &(3) 若f(b)f'(b)>0,取x_0=b,否則取x_0=a;\\ &(4) x_1=x_0-\frac{f(x_0)}{f'(x_0)};\\ &(5) 若|x_1-x_0|< ε, 則 則 輸出近似解x_1,否則令x_0=x_1並返回步驟4. \end{align*} \]
弦截法
就是用弦的斜率代替牛頓法中的
\[\begin{equation*} \begin{aligned} x_{n+1}&=x_{n}-\frac{f(x_{n})}{\frac{f(x_{n})-f(x_{n-1})}{x_n-x_{n-1}}}\\ &=\frac{x_{n-1}*f(x_{n})-x_{n}*f(x_{n-1})}{f(x_{n})-f(x_{n-1})} &其中n=1,2,3,4…… \end{aligned} \end{equation*} \]
\[\begin{align*} &弦截法的具體步驟：\\ &(1)選擇迭代初值 x_0、x_1 及輸入精度指標 ε>0 \\ &(1.5)也可以利用x_1=x_0-\frac{f(x_0)}{f'(x_0)}得出x_1;\\ &(2)計算迭代公式 x_{n+1}=\frac{x_{n-1}*f(x_{n})-x_{n}*f(x_{n-1})}{f(x_{n})-f(x_{n-1})} ,其中n=1,2,3,4… \\ &(3)如果 |x_{1}-x_{0}|< ε，輸出滿足精度的根x_{n+1}；否則 x_0=x_1,x_1=x_2並返回步驟2\\ \end{align*} \]
牛頓與弦截法的小小比較
容易看出兩者每一次迭代中所利用的直線方程分別為
fzero:基於牛頓法的找零點
x=fzero(f,x0) %求f(x)=0在點x0附近的零點
x=fzero(f,[a,b]) %求f(x)=0在[a,b]內的零點
>> clear
>> f='exp(x)-1';
>> x=fzero(f,-5)
x =
-1.2159e-17
>> x=fzero(f,3)
x =
3.7104e-17
>> x=fzero(f,0)
x =
0
fsolve:指定初始點求函數零點
[x,f,h]=fsolve(f,x0) %x為近似零點,f為該點處函數值,h輸出值大於零表示結果可靠,否則不可靠
>> [x,f,h]=fsolve('sin(x)-1',0)
Equation solved.
fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.
<stopping criteria details>
x =
1.5622
f =
-3.6930e-05
h =
1