[已解决] 如何求两个参数方程组联立后的公共参数?

[复制链接]
科研汪 发表于 2022-11-19 10:24:17
本帖最后由 科研汪 于 2022-11-19 10:25 编辑

二者的参数方程组都含有theta这个参数,现在要求出两个参数方程联立后的共同的theta
代码如下:
n1 = 1;p1 = 30.9662;k1 = 0.15;thetab0 = 0.1942;alphan = pi/9;ha = 3;
[theta,rR,Ra,thetaR,thetaa] = get_c(p1,n1,k1,thetab0,alphan,ha);
% 两个参数方程联立求解
function [theta,rR,Ra,thetaR,thetaa] = get_c(p,n,k,thetab0,alphan,ha)
[r,tau,S] = return_r_tau_S(p,n,k,thetab0);
% 非标准渐开线
% 直角坐标
xR = @(theta) r(theta) * cos(theta) - S(theta) .* cos(alphan) .* cos(theta + tau(theta) + alphan);
yR = @(theta) r(theta) * sin(theta) - S(theta) .* cos(alphan) .* sin(theta + tau(theta) + alphan);
% 极坐标
rR = @(theta) sqrt((xR(theta)).^2 + (yR(theta)).^2);
thetaR = @(theta) atan(yR(theta)./xR(theta));
% 高阶椭圆极坐标
Ra = @(theta) sqrt(r(theta).^2 + ha.^2  + 2 .* r(theta) .* ha .* sin(tau(theta)));
thetaa =  @(theta) theta - asin(ha * cos(tau(theta))./Ra(theta));
myfunR = @(theta) (Ra(theta) - rR(theta));
theta = fsolve(myfunR,thetab0);
end
% 参数方程中其他参数
function [r,tau,S] = return_r_tau_S(p,n,k,theta0)
r = @(theta) -p./(k*cos(n*theta) - 1);
tau = @(theta) -atan((p*abs(k*cos(n*theta) - 1).^2)./(abs(k*n*p*sin(n*theta)).*(k*cos(n*theta) - 1)));
L = @(t) (p.^2 ./ (k .* cos(n .* t) - 1).^2 + ...
    (k.^2 .* n.^2 .* p.^2 .* sin(n .* t).^2)./(k .* cos(n .* t) - 1).^4) .^(1 ./ 2);
S = @(x) integral(L,theta0,x);
end
运行结果如下:%rR,Ra,thetaR,thetaa作为验证语句,省略其函数句柄运行输出
theta =
    0.4321
验证:
rR(theta),Ra(theta),thetaR(theta),thetaa(theta)
结果:
ans =
   38.8421

ans =
   38.8421

ans =
    0.2305

ans =
    0.4265
通过分析,rR何Ra结果相同,但是thetaR和thetaR结果不同,按理说二者曲线只有一个交点,为什么求得的解的极径相同但角度不同呢?
二者曲线图如下:

交点为(38.1431,9.22877),换算为极坐标:极径ρ=39.243678123,角度φ=0.2373891233,也就是说,条件应满足:
Ra(theta)=rR(theta) = ρ,thetaR(theta)=thetaa(theta) = φ。
哪个地方出现问题了啊,求解答谢谢


最佳答案


TouAkira 发表于 2022-11-19 14:19:17
你想当然地误认为,两条参数曲线在交点处,各自参数取值应当相同,但实际上根本就没有这个关系,即并不存在所谓的公共参数
举个最简单易懂的反例,( x - 1 ) / 3 = ( y - 1 ) / ( -4 ) 与 ( x - 0 ) / 4 = ( y - 0 ) / 3 是两条直线
分别可以写成参数形式
L1 :  x1 = 4 * p1 + 0,  y1 = 3 * p1 + 0
L2 :  x2 = 3 * p2 + 1,  y2 = -4 * p2 + 1
交点应满足横纵坐标分别相等。
于是解线性方程组可得交点坐标为( 1.12, 0.84 ),但两个参数取值分别是 p1 = 0.28 和 p2 = 0.04
回到你这个算例,虽然两条曲线看起来都是以 theta 为参数(正如上例两直线也都可以认为以 p 为参数一样),但实际上应当理解为 theta_1 与 theta_2 分别独立。
  1. clear; clc; close all; % format longG;
  2. rng default;
  3. n = 1; p = 30.9662; k = 0.15; theta0 = 0.1942; alphan = pi/9; ha = 3;
  4. r = @( theta ) -p./( k*cos( n*theta ) - 1 );
  5. tau = @( theta ) -atan( ( p*abs( k*cos( n*theta ) - 1 ).^2 )./( abs( k*n*p*sin( n*theta ) ).*( k*cos( n*theta ) - 1 ) ) );
  6. L = @( t ) ( p.^2 ./ ( k .* cos( n .* t ) - 1 ).^2 + ...
  7.     ( k.^2 .* n.^2 .* p.^2 .* sin( n .* t ).^2 )./( k .* cos( n .* t ) - 1 ).^4 ) .^( 1 ./ 2 );
  8. S = @( x ) arrayfun( @( ii ) integral( L, theta0, x( ii ) ), [ 1 : 1 : numel( x ) ] );
  9. xR = @( theta ) r( theta ) .* cos( theta ) - S( theta ) .* cos( alphan ) .* cos( theta + tau( theta ) + alphan );
  10. yR = @( theta ) r( theta ) .* sin( theta ) - S( theta ) .* cos( alphan ) .* sin( theta + tau( theta ) + alphan );
  11. % 极坐标
  12. rR = @( theta ) sqrt( ( xR( theta ) ).^2 + ( yR( theta ) ).^2 );
  13. thetaR = @( theta ) atan2( yR( theta ),xR( theta ) );
  14. Ra = @( theta ) sqrt( r( theta ).^2 + ha.^2  + 2 .* r( theta ) .* ha .* sin( tau( theta ) ) );
  15. thetaa =  @( theta ) theta - asin( ha * cos( tau( theta ) )./Ra( theta ) );
  16. N = 20; % 采样点数目
  17. v_theta_Involute = linspace( 0, pi , N ); % 渐开线绘图幅角范围
  18. v_theta_Ellipse = linspace( 0, pi / 2, N ); % 椭圆绘图幅角范围
  19. v_R = rR( v_theta_Involute );
  20. v_theta_1 = thetaR( v_theta_Involute );
  21. v_Ra = Ra( v_theta_Ellipse );
  22. v_theta_2 = thetaa( v_theta_Ellipse );
  23. myfunR = @( Thetas ) norm( [ ( rR( Thetas( 1 ) ) * cos( thetaR( Thetas( 1 ) ) ) - Ra( Thetas( 2 ) ) * cos( thetaa( Thetas( 2 ) ) ) ), ...
  24.     ( rR( Thetas( 1 ) ) * sin( thetaR( Thetas( 1 ) ) ) - Ra( Thetas( 2 ) ) * sin( thetaa( Thetas( 2 ) ) ) ) ] );
  25. Sol_thetas = fsolve( myfunR, rand( 1, 2 ) );
  26. plot( v_R .* cos( v_theta_1 ), v_R .* sin( v_theta_1 ), 'r-.', ...
  27.     v_Ra .* cos( v_theta_2 ), v_Ra .* sin( v_theta_2 ), 'b-', ...
  28.     'linewidth', 1, 'MarkerSize', 5 );
  29. hold on;
  30. plot( rR( Sol_thetas( 1 ) ) * cos( thetaR( Sol_thetas( 1 ) ) ), rR( Sol_thetas( 1 ) ) * sin( thetaR( Sol_thetas( 1 ) ) ), 'kh', ...
  31.     Ra( Sol_thetas( 2 ) ) * cos( thetaa( Sol_thetas( 2 ) ) ), Ra( Sol_thetas( 2 ) ) * sin( thetaa( Sol_thetas( 2 ) ) ), 'gp', ...
  32.     'linewidth', 3/2, 'MarkerSize', 10 )
复制代码

交点坐标约为 ( 38.1435, 9.22925 )
回复此楼

2 条回复


TouAkira 发表于 2022-11-19 14:19:17
你想当然地误认为,两条参数曲线在交点处,各自参数取值应当相同,但实际上根本就没有这个关系,即并不存在所谓的公共参数
举个最简单易懂的反例,( x - 1 ) / 3 = ( y - 1 ) / ( -4 ) 与 ( x - 0 ) / 4 = ( y - 0 ) / 3 是两条直线
分别可以写成参数形式
L1 :  x1 = 4 * p1 + 0,  y1 = 3 * p1 + 0
L2 :  x2 = 3 * p2 + 1,  y2 = -4 * p2 + 1
交点应满足横纵坐标分别相等。
于是解线性方程组可得交点坐标为( 1.12, 0.84 ),但两个参数取值分别是 p1 = 0.28 和 p2 = 0.04
回到你这个算例,虽然两条曲线看起来都是以 theta 为参数(正如上例两直线也都可以认为以 p 为参数一样),但实际上应当理解为 theta_1 与 theta_2 分别独立。
  1. clear; clc; close all; % format longG;
  2. rng default;
  3. n = 1; p = 30.9662; k = 0.15; theta0 = 0.1942; alphan = pi/9; ha = 3;
  4. r = @( theta ) -p./( k*cos( n*theta ) - 1 );
  5. tau = @( theta ) -atan( ( p*abs( k*cos( n*theta ) - 1 ).^2 )./( abs( k*n*p*sin( n*theta ) ).*( k*cos( n*theta ) - 1 ) ) );
  6. L = @( t ) ( p.^2 ./ ( k .* cos( n .* t ) - 1 ).^2 + ...
  7.     ( k.^2 .* n.^2 .* p.^2 .* sin( n .* t ).^2 )./( k .* cos( n .* t ) - 1 ).^4 ) .^( 1 ./ 2 );
  8. S = @( x ) arrayfun( @( ii ) integral( L, theta0, x( ii ) ), [ 1 : 1 : numel( x ) ] );
  9. xR = @( theta ) r( theta ) .* cos( theta ) - S( theta ) .* cos( alphan ) .* cos( theta + tau( theta ) + alphan );
  10. yR = @( theta ) r( theta ) .* sin( theta ) - S( theta ) .* cos( alphan ) .* sin( theta + tau( theta ) + alphan );
  11. % 极坐标
  12. rR = @( theta ) sqrt( ( xR( theta ) ).^2 + ( yR( theta ) ).^2 );
  13. thetaR = @( theta ) atan2( yR( theta ),xR( theta ) );
  14. Ra = @( theta ) sqrt( r( theta ).^2 + ha.^2  + 2 .* r( theta ) .* ha .* sin( tau( theta ) ) );
  15. thetaa =  @( theta ) theta - asin( ha * cos( tau( theta ) )./Ra( theta ) );
  16. N = 20; % 采样点数目
  17. v_theta_Involute = linspace( 0, pi , N ); % 渐开线绘图幅角范围
  18. v_theta_Ellipse = linspace( 0, pi / 2, N ); % 椭圆绘图幅角范围
  19. v_R = rR( v_theta_Involute );
  20. v_theta_1 = thetaR( v_theta_Involute );
  21. v_Ra = Ra( v_theta_Ellipse );
  22. v_theta_2 = thetaa( v_theta_Ellipse );
  23. myfunR = @( Thetas ) norm( [ ( rR( Thetas( 1 ) ) * cos( thetaR( Thetas( 1 ) ) ) - Ra( Thetas( 2 ) ) * cos( thetaa( Thetas( 2 ) ) ) ), ...
  24.     ( rR( Thetas( 1 ) ) * sin( thetaR( Thetas( 1 ) ) ) - Ra( Thetas( 2 ) ) * sin( thetaa( Thetas( 2 ) ) ) ) ] );
  25. Sol_thetas = fsolve( myfunR, rand( 1, 2 ) );
  26. plot( v_R .* cos( v_theta_1 ), v_R .* sin( v_theta_1 ), 'r-.', ...
  27.     v_Ra .* cos( v_theta_2 ), v_Ra .* sin( v_theta_2 ), 'b-', ...
  28.     'linewidth', 1, 'MarkerSize', 5 );
  29. hold on;
  30. plot( rR( Sol_thetas( 1 ) ) * cos( thetaR( Sol_thetas( 1 ) ) ), rR( Sol_thetas( 1 ) ) * sin( thetaR( Sol_thetas( 1 ) ) ), 'kh', ...
  31.     Ra( Sol_thetas( 2 ) ) * cos( thetaa( Sol_thetas( 2 ) ) ), Ra( Sol_thetas( 2 ) ) * sin( thetaa( Sol_thetas( 2 ) ) ), 'gp', ...
  32.     'linewidth', 3/2, 'MarkerSize', 10 )
复制代码

交点坐标约为 ( 38.1435, 9.22925 )
回复此楼

科研汪 发表于 2022-11-19 14:32:47
TouAkira 发表于 2022-11-19 14:19
你想当然地误认为,两条参数曲线在交点处,各自参数取值应当相同,但实际上根本就没有这个关系,即并不存在 ...

我理解了,谢谢解答
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表