u
是自由参数,但其范围受第三个等式限制:
sin(v)*(sqrt(1)*cosh(u)+cos(±))=sin(±)*sinh(u)
这可以重写为:
sin(v)=-sin(alpha)*sinh(u)/(sqrt(1-y)*cosh(u)+cos(alpha))
知道这一点
abs(sin(v))<=1.
给出以下条件:
u
.
使用该
cosh(x)^2-sinh(x)^2=1
条件变为:
abs(-sin(alpha)*sqrt(cosh(u)^2-1)/(sqrt(1-y)*cosh(u)+cos(alpha))<=1.
因为cosh(x)是一个偶数函数,所以上面的表达式也是。因此,只需计算
-sin(alpha)*sqrt(cosh(u)^2-1)/(sqrt(1-y)*cosh(u)+cos(alpha)<=1.
我们想知道最大值
u
表达式适用的(
u_max
),因为我们知道
u
范围有限
[-u_max,u_max]
. 所以我们需要解决
-sin(alpha)*sqrt(cosh(u)^2-1)/(sqrt(1-y)*cosh(u)+cos(alpha)=1
这是一个二阶多项式,因此有2个解。我们对实解感兴趣,如果所有解都是虚解,那么
u
.
将其放入MATLAB,将生成以下代码:
g = [0 .05 .15 .2]; % different gammas
p = {'--k' ':g' '-r' '-b'}; % for plotting
alpha = 8*pi/5;
syms x
for i=1:length(g)
gamma = g(i);
% solve condition
sinv = -sin(alpha).*sinh(x)./(sqrt(1-gamma).*cosh(x)+cos(alpha));
sols = solve((sinv) == 1, x); % will have max 2 solutions
% pick right solution
if isreal(sols(1))
u_max = double(sols(1));
elseif isreal(sols(2))
u_max = double(sols(2));
else % both sols imaginary: no limit on u_max
u_max = 5;
end
u = -u_max:0.003:u_max;
sinv = -sin(alpha).*sinh(u)./(sqrt(1-gamma).*cosh(u)+cos(alpha));
cosv = sqrt(1-sinv.^2); % actually +-sqrt(), taken into account when plotting
px = cosh(u).*cosv;
py = sinh(u).*sinv;
plot(px,py, p{i},-px,py, p{i})
hold on
end
hold off
编辑:更新代码
g = [0 .01 .1 .75];
p = {'--k' ':g' '-r' 'ob'};
alpha = 4*pi/3;
syms x
for i=1:4
gamma = g(i);
interval = inf;
sinv = -sin(alpha).*sinh(x)./(sqrt(1-gamma).*cosh(x)+cos(alpha));
sols = solve((sinv) == 1, x); % will have max 2 solutions
if length(sols) > 1
if isreal(sols(1)) && isreal(sols(2)) % if there are 2 real solutions, interval between is valid or unvalid
if eval(subs(sinv,x,(double(sols(1))+double(sols(2)))/2)) > 1 %interval inbetween is unvalid => u ok everywhere except in interval
u_max = double(min(sols(1),sols(2)));
u_min = double(max(sols(1),sols(2)));
interval = 0;
else %interval inbetween is valid => u ok in interval
u_max = double(max(sols(1),sols(2)));
u_min = double(min(sols(1),sols(2)));
interval = 1;
end
elseif isreal(sols(1))
u_max = double(sols(1));
elseif isreal(sols(2))
u_max = double(sols(2));
else
u_max = 3;
end
elseif isreal(sols)
if eval(subs(sinv,x,sols-.1)) < 1 && eval(subs(sinv,x,sols+.1)) < 1
u_max = 3;
else
u_max = double(sols);
end
elseif eval(subs(sinv,x,1)) < 1
u_max = 3;
else
u_max = 0;
end
if interval == 1
u1 = u_min:0.003:u_max;
u2 = -u1;
sinv1 = -sin(alpha).*sinh(u1)./(sqrt(1-gamma).*cosh(u1)+cos(alpha));
sinv2 = -sin(alpha).*sinh(u2)./(sqrt(1-gamma).*cosh(u2)+cos(alpha));
cosv1 = sqrt(1-sinv1.^2);
cosv2 = sqrt(1-sinv2.^2);
if imag(cosv1) < 10^(-6)
cosv1 = real(cosv1);
end
if imag(cosv2) < 10^(-6)
cosv2 = real(cosv2);
end
if imag(cosv1) < 10^(-6)
cosv1 = real(cosv1);
end
if imag(cosv2) < 10^(-6)
cosv2 = real(cosv2);
end
px1 = cosh(u1).*cosv1;
py1 = sinh(u1).*sinv1;
px2 = cosh(u2).*cosv2;
py2 = sinh(u2).*sinv2;
plot(([-px1(end:-1:1) px1]),([py1(end:-1:1) py1]), p{i}, ([-px2(end:-1:1) px2]),([py2(end:-1:1) py2]), p{i})
hold on
elseif interval == 0
u1 = -u_max:0.003:u_max;
u2 = u_min:0.003:3;
sinv1 = -sin(alpha).*sinh(u1)./(sqrt(1-gamma).*cosh(u1)+cos(alpha));
sinv2 = -sin(alpha).*sinh(u2)./(sqrt(1-gamma).*cosh(u2)+cos(alpha));
cosv1 = sqrt(1-sinv1.^2);
cosv2 = sqrt(1-sinv2.^2);
if imag(cosv1) < 10^(-6)
cosv1 = real(cosv1);
end
if imag(cosv2) < 10^(-6);
cosv2 = real(cosv2);
end
px1 = cosh(u1).*cosv1;
py1 = sinh(u1).*sinv1;
px2 = cosh(u2).*cosv2;
py2 = sinh(u2).*sinv2;
plot([-px1(end:-1:1) (px1)],[py1(end:-1:1) (py1)], p{i}, ([-px2(end:-1:1) px2]),([py2(end:-1:1) py2]), p{i})
hold on
else
u = -u_max:0.003:u_max;
sinv1 = -sin(alpha).*sinh(u)./(sqrt(1-gamma).*cosh(u)+cos(alpha));
cosv1 = sqrt(1-sinv1.^2);
if imag(cosv1) < 10^(-6)
cosv1 = real(cosv1);
end
px1 = cosh(u).*cosv1;
py1 = sinh(u).*sinv1;
plot(-px1(end:-1:1), py1(end:-1:1), p{i}, px1, py1, p{i})
hold on
end
end
hold off
