具有两个参数和一个函数的参数图

Parametric Plot with two parameter and a function

我想在 Matlab 中使用这些方程来绘制这个参数方程

p_x = p_0*coshu*cosv, p_y = p_0*sinhu*sinv

sinv*( sqrt(1 − γ)*coshu + cosα) = −sinα *sinhu

我需要 p_y/p_0 vs p_x/p_0 之间的情节。如图

当'α = 8*pi/5时,u 是一个自由参数。和 γ = 0, 0.05, 0.15, 0.2

我尝试用代码求解上述方程式;

close all
clear all
clc
a = 8*pi/5        % 'a' as alpha
 %z=0;           % 'z'   as  gamma
z=0.15      
u = -5:0.003:5;
x = cosh(u).*sqrt(1 - (sin(a).*sin(a).*sinh(u).*sinh(u))./square((sqrt(1-z).*cosh(u) + cos(a))));
% where x =  p_x/p_0  and y = p_y/p_0
y= -1.*((sinh(u).*sinh(u).*sin(a))./(1*sqrt(1-z).*cosh(u) + cos(a)));
plot(x,y,'-k')

评论中解方程1的另一个尝试(sign(cos(v))仍然有问题):

clc; clear;
alpha=8*pi/5; gamma=0.05;
t=1;
Py0={};
for Px0=-3:.5:3
  syms u
  F=cosh(u)*sqrt(1 - ((sin(alpha)*sinh(u))/(sqrt(1-gamma)*cosh(u)+cos(alpha)))^2 )-Px0;
  u=double(solve(F));
  Py0{t}=sinh(u).*(-(sin(alpha).*sinh(u))./(sqrt(1-gamma).*cosh(u)+cos(alpha)));
  t=t+1;
  clear u
end;
Py0
% plot(-3:0.5:3,Py0)

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