Matlab求导

Matlab求导 怎么来对已经求出来的速度曲线进行求导, 我现在通过状态空间法，求解出来一个结构的速度和位移曲线，现在想通过对速度求导的方式来得到加速度，请各位好手指点，具体的达式是怎么样的，这是个动力系统，我希望能通过求导来画出加速度的曲线 中心差分法: 0.5*(diff(y(1:end-1)+diff(2:end))/dt 这样比原始数据的头尾少两个点,可以自己补一下 以试试根据数据用polyfit求出多项式拟合，然后再diff得到加速度表达式，再作图 本帖最近评分记录 , lxq +18 2007-6-22 11:28 理由:启发引导性回答 我用matlab求导后，用plot做曲线做不出来，总提示错误，不知道为什么，谁知道怎么出图啊,具体步骤说一下 你直接对y数据求导得到的数据长度比原来少了一个 plot肯定报错，x和y尺度不符 对符号函数求导之后得到符号函数 要么ezplot，要么代入数据再画图 matlab的求导命令与求导法 技术专区 2009-03-10 16:19 阅读543 评论1 字号: 大 中 小 ,(matlab命令. 建立符号变量命令sym和syms调用格式: x=sym('x')， 建立符号变量x; syms x y z ， 建立多个符号变量x，y，z; matlab求导命令diff调用格式: diff(函数) ， 求的一阶导数; diff(函数， n) ， 求的n阶导数(n是具体整数); diff(函数，变量名)， 求对的偏导数; diff(函数， 变量名，n) ，求对的n阶偏导数; matlab求雅可比矩阵命令jacobian，调用格式: jacobian([函数;函数; 函数]， [])给出矩阵: ,(导数概念. 导数是函数的变化率，几何意义是曲线在一点处的切线斜率. (,)点导数是一个极限值. 例1.设，用定义计算. 解:在某一点的导数定义为极限: 我们记，输入命令: syms h;limit((exp(0+h)-exp(0))/h，h，0) 得结果:ans=,(可知 (,)导数的几何意义是曲线的切线斜率. 例2.画出在处()的切线及若干条割线，观察割线的变化趋势. 解:在曲线上另取一点，则的方程是: .即 取，分别作出几条割线. h=[3，2，1，0.1，0.01];a=(exp(h)-1)./h;x=-1:0.1:3; plot(x，exp(x)，’r’);hold on for i=1:5; plot(h(i)，exp(h(i))，’r.’) plot(x，a(i)*x+1) end axis square 作出在处的切线 plot(x，x+1，’r’) 从图上看，随着与越来越接近，割线越来越接近曲线的割线( ,(求一元函数的导数. (,)的一阶导数. 例3.求的导数. 解:打开matlab指令窗，输入指令: dy_dx=diff(sin(x)/x). 得结果: dy_dx=cos(x)/x-sin(x)/x^2. matlab的函数名允许使用字母、空格、下划线及数字，不允许使用其他字符，在这里我们用dy_dx表示 例4.求的导数. 解: 输入命令: dy_dx=diff(log(sin(x))). 得结果: dy_dx=cos(x)/sin(x). 在matlab中，函数用log(x)表示，而log10(x)表示 例5.求的导数. 解: 输入命令:dy_dx=diff((x^2+2*x)^20). 得结果: dy_dx=20*(x^2+2*x)^19*(2*x+2). 注意输入时应为2*x. 例6.求的导数. 解: 输入命令: dy_dx=diff(x^x). 得结果: dy_dx =x^x*(log(x)+1). 利用matlab 命令diff一次可以求出若干个函数的导数. 例7.求下列函数的导数: 1. 2. 3. 4. 解: 输入命令: a=diff([sqrt(x^2- 2*x+5)，cos(x^2)+2*cos(2*x)，4^(sin(x))， log(log(x))]). 得结果: a= [1/2/(x^2-2*x+5)^(1/2)*(2*x-2)，-2*sin(x^2)*x-4*sin(2*x)， 4^sin(x)*cos(x)*log(4)， 1/x/log(x)]. dy1_dx=a(1) dy1_dx=1/2/(x^2-2*x+5)^(1/2)*(2*x-2). dy2_dx=a(2) dy2_dx=-2*sin(x^2)*x-4*sin(2*x). dy3_dx=a(3) dy3_dx=4^sin(x)*cos(x)*log(4). dy4_dx=a(4) dy4_dx=1/x/log(x). 由本例可以看出，matlab函数是对矩阵或向量进行操作的，a(i)表示向量a的第i个分量. (,)参数方 程所确定的函数的导数. 设参数方程确定函数，则的导数 例8.设，求 解: 输入命令: dx_dt=diff(a*(t-sin(t)));dy_dt=diff(a*(1-cos(t))); dy_dx=dy_dt/dx_dt. 得结果: dy_dx=sin(t)/(1-cos(t)). 其中分号的作用是不显示结果. ,(求多元函数的偏导数. 例9.设 求 u的一阶偏导数. 解: 输入命令: diff((x^2+y^2+z^2)^(1/2)， x). 得结果: ans=1/(x^2+y^2+z^2)^(1/2)*x. 在命令中将末尾的x换成y将给出y的偏导数: ans=1/(x^2+y^2+z^2)^(1/2)*y. 也可以输入命令: jacobian((x^2+y^2+z^2)^(1/2)，[x y]). 得结果: ans=[1/(x^2+y^2+z^2)^(1/2)*x， 1/(x^2+y^2+z^2)^(1/2)*y] 给出矩阵 例10.求下列函数的偏导数: 1. 2. 解: 输入命令: diff(atan(y/x). 得结果: ans=-y/x^2/(1+y^2/x^2). 输入命令: diff(atan(y/x)， y). 得结果: ans=1/x/(1+y^2/x^2). 输入命令: diff(x^y， x). 得结果: ans=x^y*y/x. 输入命令: diff(x^y， y). 得结果: ans=x^y*log(x). 使用jacobian命令求偏导数更为方便. 输入命令: jacobian([atan(y/x)，x^y]，[x，y]). 得结果: ans=[ -y/x^2/(1+y^2/x^2)， 1/x/(1+y^2/x^2)] [ x^y*y/x， x^y*log(x)]. ,(求高阶导数或高阶偏导数. 例11.设 ，求. 解:输入指令: diff(x^2*exp(2*x)，x，20). 得结果: ans = 99614720*exp(2*x)+20971520*x*exp(2*x)+1048576*x^2*exp(2*x) 例3.12.设，求,, 解:输入命令: diff(x^6-3*y^4+2*x^2*y^2，x，2) 可得到: ans=30*x^4+4*y^2. 将命令中最后一个x换为y得: ans=-36*y^2+4*x^2. 输入命令: diff(diff(x^6-3*y^4+2*x^2*y^2，x)，y) 可得: ans=8*x*y 同学们可自己计算比较它们的结果. 注意命令:diff(x^6-3*y^4+2*x^2*y^2，x，y)，是对y求偏导数，不是求 help diff If X is a vector, then diff(X) returns a vector, one element shorter than X, of differences between adjacent elements: y=f(x)，则x(i)点处的导数为(y(i+1)-y(i-1))/(x(i+1)-x(i-1)) Re:【求助】怎样用matlab计算离散曲线上各点的曲率啊? 导数用差分代替 Re:【求助】怎样用matlab计算离散曲线上各点的曲率啊? 版主说得也对，还可以用插值法先求出函数表达式，再求导 】 最优化计算之黄金搜索法matlab程序代码 Author: aiaishike 14 八 function [xo,fo]=Opt_Golden(f,a,b,TolX,TolFun,k) %%%%黄金搜索算法求在区间[a,b]上的最优化解 %f为目标函数，TolX为x阈值，TolFun为函数阈值，k为迭代次数 r =(sqrt(5)-1)/2; %r为黄金分割点值， h = b-a; %区间宽度 rh = r*h; %%%取两点c、d,并计算相应的函数值fc和fd c = b-rh; d = a+rh; fc = feval(f,c); fd = feval(f,d); %%%算法第二步判断是否停止迭代 if k <= 0 | (abs(h) < TolX & abs(fc – fd) < TolFun) if fc <= fd xo = c; fo = fc; else xo = d; fo = fd; end if k == 0 fprintf(’最好设定迭代次数大于0′); end %%%%算法第三步，进行新一轮迭代 else if fc < fd [xo,fo] = Opt_Golden(f,a,d,TolX,TolFun,k-1); else [xo,fo] = Opt_Golden(f,c,b,TolX,TolFun,k-1); end end , 0 Comments , Filed under: Matlab算法集, 专题学习资源库 偏微分方程求解有限元法之得出基函数matlab程序代码 Author: aiaishike 14 八 function Show_Basis() %N = [-1 1;1 1;1 -1;-1 -1;0.2 0.5]; %节点集合 N = [1 0;0 1;1 2;2 1;1.2 1.5]; %节点集合 N_n = size(N,1); % 总节点数 S = [1 2 5;2 3 5;3 4 5;1 4 5]; %区域集合 N_s = size(S,1); % 总区域数 figure(1), clf for s = 1:N_s nodes = [S(s,:) S(s,1)]; for i = 1:3 plot([N(nodes(i),1) N(nodes(i + 1),1)],[N(nodes(i),2) N(nodes(i+1),2)]), hold on end end text(0.8,0.6,’S1′);text(0.8,1.6,’S2′);text(1.4,1.5,’S3′),text(1.4,0.6,’S4′); %基函数 p = fem_basis_ftn(N,S); %x0 = -1; xf = 1; y0 = -1; yf = 1; %graphic region x0 = 0; xf = 2; y0 = 0; yf = 2; %graphic region figure(2), clf Mx = 50; My = 50; c = [0 1 2 3 0]; %节点的值 dx = (xf – x0)/Mx; dy = (yf – y0)/My; xi = x0 + [0:Mx]*dx; yi = y0 + [0:My]*dy; i_ns = [1 2 3 4 5]; %节点标号 for itr = 1:5 i_n = i_ns(itr); if itr == 1 for i = 1:length(xi) for j = 1:length(yi) Z(j,i) = 0; for s = 1:N_s if inpolygon(xi(i),yi(j), N(S(s,:),1),N(S(s,:),2)) > 0 Z(j,i) = p(i_n,s,1) + p(i_n,s,2)*xi(i) + p(i_n,s,3)*yi(j); break; end end end end subplot(321), mesh(xi,yi,Z);title(itr) %节点1的基函数 else c1 = zeros(size(c)); c1(i_n) = 1; subplot(320 + itr) trimesh(S,N(:,1),N(:,2),c1) %节点 2-5的基函数 title(itr); end end c = [0 1 2 3 0]; subplot(326) trimesh(S,N(:,1),N(:,2),c);title(’基函数的线型组合’) , 0 Comments , Filed under: Matlab算法集, 专题学习资源库 偏微分方程求解有限元法之为每一个节点和子区域构造基函数matlab 程序代码 Author: aiaishike 14 八 function [U,c] = fem_coef(f,g,p,c,N,S,N_i) %p(i,s,1:3): 基函数 ftn phi_i系数 %c = [ .1 1 . 0 0 .] 边界节点取1，内节点取0 %N(n,1:2) : 第n个节点的x和y坐标 %S(s,1:3) : 第s个子区域的节点#s %N_i : 内节点个数 %U(s,1:3) : 每个区域的 p1 + p2(s)x + p3(s)y 的坐标 N_n = size(N,1); % 总共节点数 N_s = size(S,1); % 总共子区域数 d=zeros(N_i,1); N_b = N_n-N_i; for i = N_b+1:N_n for n = 1:N_n for s = 1:N_s xy = (N(S(s,1),:) + N(S(s,2),:) + N(S(s,3),:))/3; %重心 p_vctr = [p([i n],s,1) p([i n],s,2) p([i n],s,3)]; tmpg(s) = sum(p(i,s,2:3).*p(n,s,2:3))-g(xy(1),xy(2))*p_vctr(1,:)*[1 xy]‘*p_vctr(2,:)*[1 xy]‘; dS(s) = det([N(S(s,1),:) 1; N(S(s,2),:) 1;N(S(s,3),:) 1])/2; %子区域 if n == 1, tmpf(s) = -f(xy(1),xy(2))*p_vctr(1,:)*[1 xy]‘; end end A12(i – N_b,n) = tmpg*abs(dS)’; end d(i-N_b) = tmpf*abs(dS)’; end d = d – A12(1:N_i,1:N_b)*c(1:N_b)’; c(N_b + 1:N_n) = A12(1:N_i,N_b+1:N_n)\d; for s = 1:N_s for j = 1:3, U(s,j) = c*p(:,s,j); end end , 0 Comments , Filed under: Matlab算法集, 专题学习资源库 偏微分方程求解有限元法之产生基函数matlab程序代码 Author: aiaishike 14 八 function p = fem_basis_ftn(N,S) %p(i,s,1:3): 基函数 ftn phi_i系数 % 共s个子区域 %N(n,1:2) : 第n个节点的x和y坐标 %S(s,1:3) : 第s个子区域的节点#s N_n = size(N,1); %总共节点数 N_s = size(S,1); % 总子区域个数 for n = 1:N_n for s = 1:N_s for i = 1:3 A(i,1:3) = [1 N(S(s,i),1:2)]; b(i) = (S(s,i) == n); %The nth basis ftn is 1 only at node n. end pnt=A\b’; for i=1:3, p(n,s,i) = pnt(i); end end end , 0 Comments , Filed under: Matlab算法集, 专题学习资源库 显示中心差分法求解二维双曲线方程matlab程序代码 Author: aiaishike 14 八 function [u,x,y,t] = Wave2(A,D,T,it0,i1t0,bxyt,Mx,My,N) %解方程 a(u_xx + u_yy) = u_tt for D(1) <= x <= D(2), D(3) <= y <= D(4), 0 <= t <= T % 初始条件: u(x,y,0) = it0(x,y), u_t(x,y,0) = i1t0(x,y) % 边界条件: u(x,y,t) = bxyt(x,y,t) for (x,y) on Boundary % Mx/My x轴和y轴的等分段数 % N = 时间轴的等分段数 dx = (D(2)-D(1))/Mx; x = D(1)+[0:Mx]*dx; dy = (D(4)- D(3))/My; y = D(3)+[0:My]‘*dy; dt = T/N; t = [0:N]*dt; %初始化 u = zeros(My+1,Mx + 1); ut = zeros(My + 1,Mx + 1); for j = 2:Mx for i = 2:My u(i,j) = it0(x(j),y(i)); ut(i,j) = i1t0(x(j),y(i)); end end adt2 = A*dt*dt; rx = adt2/(dx*dx); ry = adt2/(dy*dy); rxy1 = 1- rx – ry; rxy2 = rxy1*2; u_1 = u; for k = 0:N t = k*dt; for i = 1:My + 1 %边界条件 u(i,[1 Mx + 1]) = [bxyt(x(1),y(i),t) bxyt(x(Mx + 1),y(i),t)]; end for j = 1:Mx + 1 u([1 My + 1],j) = [bxyt(x(j),y(1),t); bxyt(x(j),y(My + 1),t)]; end if k == 0 for i = 2:My for j = 2:Mx %(11.3.8) u(i,j) = 0.5*(rx*(u_1(i,j – 1) + u_1(i,j + 1))… + ry*(u_1(i – 1,j)+u_1(i + 1,j))) + rxy1*u(i,j) + dt*ut(i,j); end end else for i = 2:My for j = 2:Mx %11.3.7 u(i,j) = rx*(u_1(i,j – 1)+ u_1(i,j + 1))+ ry*(u_1(i – 1,j) + u_1(i + 1,j)) + rxy2*u(i,j) -u_2(i,j); end end end u_2 = u_1; u_1 = u; end , 0 Comments , Filed under: Matlab算法集, 专题学习资源库 双曲线偏微分方程求解之显示中心差分法matlab程序代码 Author: aiaishike 14 八 function [u,x,t] = ECD_Wave(A,xf,T,it0,i1t0,bx0,bxf,M,N) %解方程a u_xx = u_tt for 0<=x<=xf, 0<=t<=T % 初始条件: u(x,0) = it0(x), u_t(x,0) = i1t0(x) % 边界条件: u(0,t)= bx0(t), u(xf,t) = bxf(t) % M :沿x轴的等分段数 % N :沿y轴的等分段数 dx = xf/M; x = [0:M]‘*dx; dt = T/N; t = [0:N]*dt; for i = 1:M + 1 u(i,1) = it0(x(i)); end for k = 1:N + 1 u([1 M + 1],k) = [bx0(t(k)); bxf(t(k))]; end r = A*(dt/dx)^ 2; r1 = r/2; r2 = 2*(1 – r); u(2:M,2) = r1*u(1:M – 1,1) + (1 – r)*u(2:M,1) + r1*u(3:M + 1,1) + dt*i1t0(x(2:M)); %(11.3.4) for k = 3:N + 1 u(2:M,k) = r*u(1:M – 1,k – 1) + r2*u(2:M,k-1) + r*u(3:M + 1,k – 1)- u(2:M,k – 2); %(11.3.3) end