学习matlab之第三周
学习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
