【matlab基础知识代码】（七）单变量函数求导&偏导数的计算

最新推荐文章于 2024-07-19 21:51:00 发布

庄宝007

最新推荐文章于 2024-07-19 21:51:00 发布

阅读量967

点赞数 16

文章标签： matlab 开发语言

本文链接：https://blog.csdn.net/2301_80118808/article/details/138442468

版权

>> syms x; f=sin(x)/(x^2+4*x+3); f1=diff(f)
 
f1 =
 
cos(x)/(x^2 + 4*x + 3) - (sin(x)*(2*x + 4))/(x^2 + 4*x + 3)^2

>>  ezplot(f,[0,5]), hold on; ezplot(f1,[0,5])
>> f2=diff(f,4)
 
f2 =
 
sin(x)/(x^2 + 4*x + 3) + (12*sin(x))/(x^2 + 4*x + 3)^2 + (24*sin(x))/(x^2 + 4*x + 3)^3 - (24*cos(x)*(2*x + 4)^3)/(x^2 + 4*x + 3)^4 - (12*sin(x)*(2*x + 4)^2)/(x^2 + 4*x + 3)^3 - (48*sin(x)*(2*x + 4)^2)/(x^2 + 4*x + 3)^4 + (24*sin(x)*(2*x + 4)^4)/(x^2 + 4*x + 3)^5 + (4*cos(x)*(2*x + 4))/(x^2 + 4*x + 3)^2 + (16*cos(x)*(2*x + 4))/(x^2 + 4*x + 3)^3 + (8*cos(x)*(8*x + 16))/(x^2 + 4*x + 3)^3 - (6*sin(x)*(2*x + 4)*(8*x + 16))/(x^2 + 4*x + 3)^4

蓝色的函数原函数，红色的线导数曲线

>>  tic, diff(f,x,50); toc
历时 1.071734 秒。

>> syms t f(t); G=simplify(diff(t^2*sin(t)*f,t,3))
 
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)
 
>> simplify(subs(G,f,exp(-t))), simplify(diff(t^2*sin(t)*exp(-t),3)-ans)
 
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

>> syms t;f=t^2*sin(t)*exp(-t);diff(f,3)

ans =

6*exp(-t)*cos(t) - 6*exp(-t)*sin(t) + 2*t^2*exp(-t)*cos(t) + 2*t^2*exp(-t)*sin(t) - 12*t*exp(-t)*cos(t)

>> 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)
 
H =
 
[     4*sin(5*x),               exp(-4*x^2)]
[3*x^2 + 4*x + 1, 2^(1/2)*(2*x^2 + 1)^(1/2)]
 
 
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)]

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

>> 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

偏导数的计算

>> 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)

>> [x0,y0]=meshgrid(-3:.2:2,-2:.2:2); z0=double(subs(z,{x,y},{x0,y0})); surf(x0,y0,z0), zlim([-0.7 1.5])

1. `meshgrid(-3:.2:2,-2:.2:2)`：这行代码创建了一个二维网格，其中 x 取值从 -3 到 2，y 取值从 -2 到 2，步长为 0.2。这将生成两个矩阵 x0 和 y0，分别包含了 x 和 y 的所有组合。

2. `z0=double(subs(z,{x,y},{x0,y0}))`：这行代码计算了一个函数 z 在给定 x 和 y 值的情况下的值。这里假设之前已经定义了一个符号表达式 z，其中包含了 x 和 y。通过 subs 函数，用 x0 和 y0 替换 z 中的 x 和 y，并将结果转换为双精度数值。

3. `surf(x0,y0,z0), zlim([-0.7 1.5])`：这行代码使用 surf 函数绘制了一个三维曲面图，其中 x0 和 y0 为网格坐标，z0 为对应的函数值。然后使用 zlim 函数设置了 z 轴的显示范围为 -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绘制等高线图并在图上添加矢量场的功能。首先，contour(x0, y0, z0, 30)用于绘制二维等高线图，其中x0和y0是网格的x和y坐标，z0是对应的高度值，30表示绘制30个等高线。接着，hold on用于保持当前图形以便在其上添加其他图形。然后，zx0和zy0分别计算了在点(x0, y0)处的梯度值。其中，subs(zx, {x, y}, {x0, y0})和subs(zy, {x, y}, {x0, y0})用于计算梯度的x和y分量。最后，quiver(x0, y0, -zx0, -zy0)用于在等高线图上绘制矢量场，其中-xz0和-zy0是梯度的反向，表示矢量的方向。

>> syms x y z; f=sin(x^2*y)*exp(-x^2*y-z^2); df=diff(diff(diff(f,x,2),y),z);
df=simplify(df)
 
df =
 
-4*z*exp(- y*x^2 - z^2)*(cos(x^2*y) - sin(x^2*y) + 4*x^4*y^2*cos(x^2*y) + 4*x^4*y^2*sin(x^2*y) - 10*x^2*y*cos(x^2*y))

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

这个函数是用来计算一个函数在某一点处的n阶偏导数的差分逼近值的。参数f是要计算的函数，x和y是自变量，n是要计算的阶数。首先检查n是否为正整数，如果不是则报错。然后计算函数f对x和y的偏导数的比值F1，并将其赋给dy。然后通过循环计算n阶偏导数的逼近值，每次迭代都使用之前计算得到的F1和上一阶的偏导数值。最后返回计算得到的n阶偏导数的逼近值dy。

>>  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))
 
>>  F2=impldiff(f,x,y,2), F3=impldiff(f,x,y,3), [n,d]=numden(F3), simplify(n)
 
F2 =
 
-(2*(3*x^6 + 3*x^5*y - 12*x^5 + 3*x^4*y^2 - 12*x^4*y + 7*x^4 - 12*x^3*y^2 + 16*x^3*y + 16*x^3 + 16*x^2*y^2 - 8*x^2*y - 6*x^2 - 8*x*y^2 + 8*x*y - 8*x + 8*y^2 + 4))/(x^2*(x + 2*y)^3*(x - 2)^2)
 
 
F3 =
 
(2*(- 27*x^9 - 27*x^8*y + 162*x^8 - 27*x^7*y^2 + 162*x^7*y - 241*x^7 + 162*x^6*y^2 - 308*x^6*y - 204*x^6 - 276*x^5*y^2 + 132*x^5*y + 582*x^5 + 64*x^4*y^3 + 36*x^4*y^2 + 216*x^4*y + 64*x^4 + 32*x^3*y^4 - 192*x^3*y^3 + 408*x^3*y^2 - 208*x^3*y - 444*x^3 - 96*x^2*y^4 + 384*x^2*y^3 - 336*x^2*y^2 + 192*x^2*y + 48*x^2 + 192*x*y^4 - 256*x*y^3 + 192*x*y^2 - 96*x*y + 144*x - 128*y^4 - 96*y^2 - 48))/(x^3*(x + 2*y)^5*(x - 2)^3)
 
 
n =
 
- 54*x^9 - 54*x^8*y + 324*x^8 - 54*x^7*y^2 + 324*x^7*y - 482*x^7 + 324*x^6*y^2 - 616*x^6*y - 408*x^6 - 552*x^5*y^2 + 264*x^5*y + 1164*x^5 + 128*x^4*y^3 + 72*x^4*y^2 + 432*x^4*y + 128*x^4 + 64*x^3*y^4 - 384*x^3*y^3 + 816*x^3*y^2 - 416*x^3*y - 888*x^3 - 192*x^2*y^4 + 768*x^2*y^3 - 672*x^2*y^2 + 384*x^2*y + 96*x^2 + 384*x*y^4 - 512*x*y^3 + 384*x*y^2 - 192*x*y + 288*x - 256*y^4 - 192*y^2 - 96
 
 
d =
 
x^3*(x + 2*y)^5*(x - 2)^3
 
 
ans =
 
- 54*x^9 - 54*x^8*y + 324*x^8 - 54*x^7*y^2 + 324*x^7*y - 482*x^7 + 324*x^6*y^2 - 616*x^6*y - 408*x^6 - 552*x^5*y^2 + 264*x^5*y + 1164*x^5 + 128*x^4*y^3 + 72*x^4*y^2 + 432*x^4*y + 128*x^4 + 64*x^3*y^4 - 384*x^3*y^3 + 816*x^3*y^2 - 416*x^3*y - 888*x^3 - 192*x^2*y^4 + 768*x^2*y^3 - 672*x^2*y^2 + 384*x^2*y + 96*x^2 + 384*x*y^4 - 512*x*y^3 + 384*x*y^2 - 192*x*y + 288*x - 256*y^4 - 192*y^2 - 96

>> syms x y;f=x^2+x*y+y^2-3;
F1=impldiff(f,x,y,1), f2=impldiff(f,x,y,2); 
F2=subs(f2,x^2+x*y+y^2,3), f3=impldiff(f,x,y,3); 
F3=subs(f3,x^2+x*y+y^2,3), f4=impldiff(f,x,y,4); 
F4=subs(f4,x^2+x*y+y^2,3) 
 
F1 =
 
-(2*x + y)/(x + 2*y)
 
 
F2 =
 
-18/(x + 2*y)^3
 
 
F3 =
 
-(162*x)/(x + 2*y)^5
 
 
F4 =
 
-(216*(4*x^4 + 5*x^3*y + 6*x^2*y^2 + 2*x*y^3 + y^4))/(x + 2*y)^7

>>  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]

>> syms x y; f=(x^2-2*x)*exp(-x^2-y^2-x*y); H=simplify(hessian(f,[x,y]))

H =

[exp(- x^2 - x*y - y^2)*(4*x^4 + 4*x^3*y - 8*x^3 + x^2*y^2 - 8*x^2*y - 10*x^2 - 2*x*y^2 - 4*x*y + 12*x + 4*y + 2), exp(- x^2 - x*y - y^2)*(2*x^4 + 5*x^3*y - 4*x^3 + 2*x^2*y^2 - 10*x^2*y - 3*x^2 - 4*x*y^2 - 4*x*y + 4*x + 4*y)]
[ exp(- x^2 - x*y - y^2)*(2*x^4 + 5*x^3*y - 4*x^3 + 2*x^2*y^2 - 10*x^2*y - 3*x^2 - 4*x*y^2 - 4*x*y + 4*x + 4*y), x*exp(- x^2 - x*y - y^2)*(x - 2)*(x^2 + 4*x*y + 4*y^2 - 2)]

早期方法

>> X=[x,y];H=jacobian(jacobian(f,X),X)
 
H =
 
[2*exp(- x^2 - x*y - y^2) + 2*exp(- x^2 - x*y - y^2)*(- x^2 + 2*x) - exp(- x^2 - x*y - y^2)*(- x^2 + 2*x)*(2*x + y)^2 - 2*exp(- x^2 - x*y - y^2)*(2*x - 2)*(2*x + y), exp(- x^2 - x*y - y^2)*(- x^2 + 2*x) - exp(- x^2 - x*y - y^2)*(2*x - 2)*(x + 2*y) - exp(- x^2 - x*y - y^2)*(- x^2 + 2*x)*(x + 2*y)*(2*x + y)]
[                       exp(- x^2 - x*y - y^2)*(- x^2 + 2*x) - exp(- x^2 - x*y - y^2)*(2*x - 2)*(x + 2*y) - exp(- x^2 - x*y - y^2)*(- x^2 + 2*x)*(x + 2*y)*(2*x + y),                                                    2*exp(- x^2 - x*y - y^2)*(- x^2 + 2*x) - exp(- x^2 - x*y - y^2)*(- x^2 + 2*x)*(x + 2*y)^2]