查看: 177|回复: 6|关注: 0

[已解决] MATLAB分段函数求偏导数,报错除数为0

[复制链接]

新手

9 麦片

财富积分


050


4

主题

21

帖子

0

最佳答案
急求:用MATLAB函数求分段函数的偏导数,运行后,问题出现在Ab(i)=subs(dAb3,e,X);处,报错除数为0,Z中没有值,求高人指点!
function diffAbc=Abc()
clear all;
clc;

%定义自变量e,n
syms e n;

%初值
l=3.4;
r=0.2;
r1=0.5;
r2=0.2;
D=11;
h=0.8;
L=26;
E=0.8;
st=1.0;
B=pi./n;

%定义e的范围
e1=l*sin(E*B)-r;
e2=D/2-l-r;
e3=sqrt(l^2+D^2/4-l*D*cos(E*B))-r;

%定义计算中n和e的取值范围,重新定义B
n=4:1:6;
e0=0:0.01:2;
[X,Y]=meshgrid(e0,n);
B=pi./Y;




%三阶段面积计算公式
Ab1=2.*Y.*l.*L.*((((sin(E.*B))/(sin(st./2)))+(l-E).*B+h./l+((r+e)./l).*(pi./2+B-tan(st./4)-csc(st./2)))-(2.*((r2-e)./l).*(tan(st./4)-st./4)-((r1-e)./l).*(tan(st./2)-st./2)));
Ab2=2.*Y.*l.*L.*(l-E).*B+(r+e).*(B+asin((l.*sin(E.*B))./(r+e)));
Ab3=2.*Y.*L.*(r+e).*(acos((l.^2+(r+e).^2-(D/2).^2)/(2.*l.*(r+e)))-pi+asin((l.*sin(E.*B))/(r+e))+E.*B);

%求偏导
dAb1=diff(Ab1,e)
dAb2=diff(Ab2,e)
dAb3=diff(Ab3,e)

%重新定义e的范围
e1=l.*sin(E.*B)-r;
e2=D/2-l-r;
e3=sqrt(l.^2+D.^2/4-l.*D.*cos(E.*B))-r;

i=1;
Ab=zeros(3,201);

while X<=e3
    if(X<e1)
        Ab(i)=subs(dAb1,e,X);
    elseif(X<e2)
        Ab(i)=subs(dAb2,e,X);
     elseif(X<e3)
        Ab(i)=subs(dAb3,e,X);
    end
    X=X+1;i=i+1;
end

Z=Ab;
surf(X,Y,Z);
shading faceted;

论坛优秀回答者

0

主题

664

帖子

121

最佳答案
  • 关注者: 24
发表于 2020-2-10 22:34:36 | 显示全部楼层 |此回复为最佳答案
仅供参考

clear all;
clc;

%定义自变量e,n
syms e n;

%初值
l=3.4;
r=0.2;
r1=0.5;
r2=0.2;
D=11;
h=0.8;
L=26;
E=0.8;
st=1.0;
B=pi./n;

%定义e的范围
e1=l*sin(E*B)-r;
e2=D/2-l-r;
e3=sqrt(l^2+D^2/4-l*D*cos(E*B))-r;

%定义计算中n和e的取值范围,重新定义B
n=4:1:6;
e0=0.01:0.01:2;
[X,Y]=meshgrid(e0,n);
B=pi./Y;




%三阶段面积计算公式
Ab1=2.*Y.*l.*L.*((((sin(E.*B))/(sin(st./2)))+(l-E).*B+h./l+((r+e)./l).*(pi./2+B-tan(st./4)-csc(st./2)))-(2.*((r2-e)./l).*(tan(st./4)-st./4)-((r1-e)./l).*(tan(st./2)-st./2)));
Ab2=2.*Y.*l.*L.*(l-E).*B+(r+e).*(B+asin((l.*sin(E.*B))./(r+e)));
Ab3=2.*Y.*L.*(r+e).*(acos((l.^2+(r+e).^2-(D/2).^2)/(2.*l.*(r+e)))-pi+asin((l.*sin(E.*B))/(r+e))+E.*B);

%求偏导
dAb1=diff(Ab1,e)
dAb2=diff(Ab2,e)
dAb3=diff(Ab3,e)

%重新定义e的范围
e1=l.*sin(E.*B)-r;
e2=D/2-l-r;
e3=sqrt(l.^2+D.^2/4-l.*D.*cos(E.*B))-r;

i=1;
Ab=zeros(9,201);

while X(:,i)<=e3(:,i)
    if(X(:,i)<e1(:,i))
        Ab(:,i)=subs(dAb1(:,i),e,X(:,i));
    elseif(X<e2)
        Ab(:,i)=subs(dAb2(:,i),e,X(:,i));
     elseif(X<e3(i))
        Ab(:,i)=subs(dAb3(:,i),e,X(:,i));
    end

    X=X+1;i=i+1;
end

Z=Ab(1:3:7,1:200);
surf(X,Y,Z);
shading faceted;
微信图片_20200210223335.png

新手

9 麦片

财富积分


050


4

主题

21

帖子

0

最佳答案
 楼主| 发表于 2020-2-11 11:26:56 | 显示全部楼层
20141303 发表于 2020-2-10 22:34
仅供参考

clear all;

Z=Ab(1:3:7,1:200)能帮忙解释一下这句是什么意思吗?

论坛优秀回答者

0

主题

664

帖子

121

最佳答案
  • 关注者: 24
发表于 2020-2-11 12:24:54 | 显示全部楼层
Z矩阵是由Ab矩阵的1,4,7行和1到200列组成,即Z为3*200矩阵

新手

9 麦片

财富积分


050


4

主题

21

帖子

0

最佳答案
 楼主| 发表于 2020-2-13 12:55:27 | 显示全部楼层
20141303 发表于 2020-2-11 12:24
Z矩阵是由Ab矩阵的1,4,7行和1到200列组成,即Z为3*200矩阵

这里面的行个列是怎么选出来的?

论坛优秀回答者

0

主题

664

帖子

121

最佳答案
  • 关注者: 24
发表于 2020-2-13 15:05:36 | 显示全部楼层
随意选的,具体根据的你的需要调整

新手

9 麦片

财富积分


050


4

主题

21

帖子

0

最佳答案
 楼主| 发表于 2020-2-16 13:01:30 | 显示全部楼层
20141303 发表于 2020-2-13 15:05
随意选的,具体根据的你的需要调整

OK,明白了,谢谢
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /4 下一条

快速回复 返回顶部 返回列表