大家都提供了道,我来提供术。我最近在自学python,题主如果想学python,可以私信微信号给我,我们一起努力啊。
我采用的是 @vczh 的算法,期间程序有问题,我手算了一条内公切线的斜率,用的是 @程德华 的方法。(我有多想念老师们出题的数值!老师们节日快乐。)
先说结论:python可以算,但是很麻烦,而且我的程序最终还是有问题。推荐使用matlab。
——————————分割线——————————————
点到直线的距离:
(还复习了Latex,真是怀念大学时光)
把分数上下除以B就可以得到k,b的等式。
圆以圆心加半径的方式传入:
圆心到直线距离等于半径,我们可以得到两个等式。
我们得到了一个二元二次方程组,最多的时候会有四组解。
在matlab里有solve函数,可以直接求出多元多次方程组的解,这样整个过程就结束了。
————————matlab党止步————————————————
python中有科学计算三驾马车,numpy,scipy和maplotlib。其中scipy的optimize中有fsolve函数,可以解多元多次方程组,但是,他需要给出初值,并且返回该初值附近的一组解。
是的,只有一组解。我所想到的办法就是从头开始赋初值,尽量别错过每一种可能性。这样程序就变得臃肿,而且效率低下。fsolve本来就是计算量巨大的函数,还要计算那么多次。(我多想有人过来给我一巴掌,笨蛋,**模块都没听过?)
我既然是来提供术的,我会把详细的过程写下来。
核心模块,解二元二次方程组。
其中[1,1]是初值,这样就可以返回[1,1]附近的方程组的解。当然我们得要四组解(如果有的话)
如何去除重复解?
PRICISION是个常数,我定义成了0.01。
可以想象整个程序的复杂度。但是根据我的观察,发现有几个特点可以利用:
1.四条公切线的斜率肯定不同,因此判定时只需要判定k即可。(半径不同的时候)
2.在其中一个循环中找到一组k,b,就可以结束内循环,跳到下一个k。
3.用很小的初值找到一组解,可以认为其他的解比刚求的解大。
4.当取得四组解时,整个过程就可以结束了。
注:1.斜率和截距超出范围外就无能为力了。
2.如果两个圆是分离的,提高斜率和截距的上下限不会影响太多性能。
其中try是为了得到第一组值的时候,result是空列表,遍历会出错。
定义一个函数funcForSolve,输入两个圆,以列表的形式返回方程组即可。
——————————回答完毕————————————————
我感觉题主应该是python新手,以下是更详细的过程。
思考一个程序的结构,我的习惯是从核心模块开始,往上和往下去思考。
——————
往上:
核心模块需要什么?一个方程组表达式的列表
如何得到这个列表?把圆的坐标带入公式
要是有圆内含怎么办?需要根据圆心距离判断
需要一个计算点到点距离的模块。
往下:
核心模块输出了什么?一个列表,里面k,b组成的tuple
用k,b组成直线的方程式输出
绘图。
绘图需要什么?
圆的坐标,直线的方程,坐标范围。
————————
我喜欢混沌的思维,总是尽可能把程序砸碎,看看有没有可以重用的模块。根据k,b返回直线方程就是可以重用的模块。
python中不支持表达式,只能把要计算的等式当作字符串存储。circle_to_func就是把圆的坐标带入公式中,而funcForSolve则是要把字符串转换成方程组。而给fsolve的参数函数最好是只有一个输入,方便输入初值。我把内含圆的判断也放到了这里。
然后就是得到k,b后,转换成方程式输出。
注意这里的方程式是要用来计算的,因此表达式得是符合python运算法则的。
然后画圆
plt是pyplot的简称,是matplotlib中的模块
画线
我根据圆心和半径建立了一个函数,返回坐标的范围。也就是上面的axes_limits()
然后就是简单的设置,显示或者保存图像。
当输入是:
输出是
以及:
我采用的是 @vczh 的算法,期间程序有问题,我手算了一条内公切线的斜率,用的是 @程德华 的方法。(我有多想念老师们出题的数值!老师们节日快乐。)
先说结论:python可以算,但是很麻烦,而且我的程序最终还是有问题。推荐使用matlab。
——————————分割线——————————————
点到直线的距离:
把分数上下除以B就可以得到k,b的等式。
![d=\frac{|kx_0-y_0+b|}{\sqrt{k^2+1}}](https://i-blog.csdnimg.cn/blog_migrate/ad24269ca43415f9db000839bd5e68fc.png)
圆以圆心加半径的方式传入:
圆心到直线距离等于半径,我们可以得到两个等式。
![\frac{|kx_1-y_1+b|}{\sqrt{k^2+1}}=r_1\\\frac{|kx_2-y_2+b|}{\sqrt{k^2+1}}=r_2](https://i-blog.csdnimg.cn/blog_migrate/832bbdf01ed528d9eb56a0a83405fda2.png)
我们得到了一个二元二次方程组,最多的时候会有四组解。
在matlab里有solve函数,可以直接求出多元多次方程组的解,这样整个过程就结束了。
————————matlab党止步————————————————
python中有科学计算三驾马车,numpy,scipy和maplotlib。其中scipy的optimize中有fsolve函数,可以解多元多次方程组,但是,他需要给出初值,并且返回该初值附近的一组解。
是的,只有一组解。我所想到的办法就是从头开始赋初值,尽量别错过每一种可能性。这样程序就变得臃肿,而且效率低下。fsolve本来就是计算量巨大的函数,还要计算那么多次。(我多想有人过来给我一巴掌,笨蛋,**模块都没听过?)
我既然是来提供术的,我会把详细的过程写下来。
核心模块,解二元二次方程组。
from scipy.optimize import fsolve
k,b = fsolve(funcForSolve, [1,1])
from numpy import *
klist = linspace(-50, 50, 200)
blist = linspace(-100, 100, 500)
for k in klist:
for b in blist:
(solved_k, solved_b) = fsolve(funcForSolve, [k,b])
if abs(k2-k1) > PRECISION and abs(b1-b2) > PRECISION:
result.append((solved_k, solved_b))
可以想象整个程序的复杂度。但是根据我的观察,发现有几个特点可以利用:
1.四条公切线的斜率肯定不同,因此判定时只需要判定k即可。(半径不同的时候)
2.在其中一个循环中找到一组k,b,就可以结束内循环,跳到下一个k。
3.用很小的初值找到一组解,可以认为其他的解比刚求的解大。
4.当取得四组解时,整个过程就可以结束了。
def solve_func():
result = []
klist = linspace(-50, 50, 200)
blist = linspace(-100, 100, 500)
(solved_k, solved_b) = (-float('Inf'), -float('Inf'))
for k in klist:
for b in blist:
if k > solved_k:
(solved_k, solved_b) = fsolve(funcForSolve, [k,b])
try:
delta_k = min([abs(solved_k - item[0]) \
for item in result])
except(IndexError, ValueError):
delta_k = 10 # for first solves
if delta_k > PRECISION:
result.append((solved_k, solved_b))
break
if len(result) >= 4:
break
return result
2.如果两个圆是分离的,提高斜率和截距的上下限不会影响太多性能。
其中try是为了得到第一组值的时候,result是空列表,遍历会出错。
定义一个函数funcForSolve,输入两个圆,以列表的形式返回方程组即可。
——————————回答完毕————————————————
我感觉题主应该是python新手,以下是更详细的过程。
思考一个程序的结构,我的习惯是从核心模块开始,往上和往下去思考。
——————
往上:
核心模块需要什么?一个方程组表达式的列表
如何得到这个列表?把圆的坐标带入公式
要是有圆内含怎么办?需要根据圆心距离判断
需要一个计算点到点距离的模块。
往下:
核心模块输出了什么?一个列表,里面k,b组成的tuple
用k,b组成直线的方程式输出
绘图。
绘图需要什么?
圆的坐标,直线的方程,坐标范围。
————————
我喜欢混沌的思维,总是尽可能把程序砸碎,看看有没有可以重用的模块。根据k,b返回直线方程就是可以重用的模块。
其中黑色箭头是程序引用的关系,红色是主函数中的顺序。我一直试图把输入控制在主函数中,后来发现函数的调用变得难以理解。
首先是怎么得到方程组。def circle_to_func(circle_1, circle_2):
(x1, y1) = circle_1[0]
r1 = circle_1[1]
(x2, y2) = circle_2[0]
r2 = circle_2[1]
dist = point_to_point(circle_1[0], circle_2[0])
assert dist > circle_1[1] and dist > circle_2[1], \
"Containing circle has no common tangent."
return '(k*%f - %f + b) ** 2 - %f**2 * (k**2 + 1)' %(x1,y1,r1),\
'(k*%f - %f + b) ** 2 - %f**2 * (k**2 + 1)' %(x2,y2,r2)
def funcForSolve(x):
'convert the string of functions to the polynomial'
k,b = x.tolist()
func = circle_to_func(CIRCLE_1, CIRCLE_2)
return [eval(func[0]), eval(func[1])]
然后就是得到k,b后,转换成方程式输出。
def argu_to_line(arguments):
result = []
for (k,b) in arguments:
result.append("%.4f*x%+.4f" %(k,b))
return result
然后画圆
fig = plt.figure(dpi=120)
ax = fig.add_subplot(111)
circle_plot_1 = plt.Circle(circle_1[0], circle_1[1], fill=False)
circle_plot_2 = plt.Circle(circle_2[0], circle_2[1], fill=False)
ax.add_artist(circle_plot_1)
ax.add_artist(circle_plot_2)
画线
limit = axes_limits(circle_1, circle_2)
x = linspace(limit[0][0], limit[0][1], 100)
for func in line_func:
plt.plot(x, eval(func), 'r')
然后就是简单的设置,显示或者保存图像。
ax.set_xlim(limit[0])
ax.set_ylim(limit[1])
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.show()
fig.savefig('common_tangent.png', transparent=True, dpi=300)
PRECISION = 0.01
CIRCLE_1 = ((10,10),1)
CIRCLE_2 = ((5,6),2)
y=0.1893*x+7.0888
y=1.0969*x-2.4528
y=0.5698*x+5.4528
y=2.3107*x-10.5888
图比较简陋,实在不好意思呐。
源代码:gt11799/Common_tangent · GitHub
原载于知乎
——————————————
github主页:https://github.com/gt11799
E-mail:gting405@163.com