根据空间中不共面的四个点坐标，求构成任意四面体的内外球

海伦公式：

四面体体积公式

六条边分别为a,b,c,a1,b1,c1。
a，b，c，a1，b1，c1，其中a与a1，b与b1，c与c1互为对边，那么有三棱锥（四面体）的体积公式为:
V=1/12sqrt(4a2b2c2-a2D2-b2E2c2F2+DEF),sqrt表示开根号，a2表示平方
D=b2+c2-a12,E=a2+c2-b12, F=a2+b2-c1**2

内切球半径公式

R = (3*V)/S s为四面体面积

最小外接球

四个点到球心的距离相等
代码如下：

import math
import numpy as np
from scipy.linalg import solve
def nball(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4):
    try:
        #六条边分别为a,b,c,a1,b1,c1
        a = math.sqrt((x1 - x2) ** 2 + (y1 - y2) ** 2 + (z1 - z2) ** 2)
        b = math.sqrt((x1 - x3) ** 2 + (y1 - y3) ** 2 + (z1 - z3) ** 2)
        c = math.sqrt((x1 - x4) ** 2 + (y1 - y4) ** 2 + (z1 - z4) ** 2)
        a1 = math.sqrt((x3 - x4) ** 2 + (y3 - y4) ** 2 + (z3 - z4) ** 2)#a的对边
        b1 = math.sqrt((x2 - x4) ** 2 + (y2 - y4) ** 2 + (z2 - z4) ** 2)#b的对边
        c1 = math.sqrt((x2 - x3) ** 2 + (y2 - y3) ** 2 + (z2 - z3) ** 2)#c的对边
        print('六条边长',[a,b,c,a1,b1,c1])
        #根据海伦公式求四面体各个面的面积，分别为s1,s2,s3,s4
        #a1b1c1边构成s1,abc1边构成s2,bca1边构成s3,acb1边构成s4
        p1 = round((a1 + b1 + c1) / 2, 3)  # 海伦公式的p值
        s1 = round(math.sqrt(p1 * (p1 - a1) * (p1 - b1) * (p1 - c1)), 3)  # 面积
        p2 = round((a + b + c1) / 2, 3)
        s2 = round(math.sqrt(p2 * (p2 - a) * (p2 - b) * (p2 - c1)), 3)
        p3 = round((b+c+a1) / 2, 3)
        s3 = round(math.sqrt(p3 * (p3 - a1) * (p3 - b) * (p3 - c)), 3)
        p4 = round((a + b1 + c) / 2, 3)
        s4 = round(math.sqrt(p4 * (p4 - a) * (p4 - b1) * (p4 - c)), 3)
        S=s1+s2+s3+s4#四面体表面积
        print('四面体表面积',round(S,3))
        #a，b，c，a1，b1，c1，其中a与a1，b与b1，c与c1互为对边，那么有三棱锥（四面体）的体积公式为:
        #V=1/12sqrt(4a**2b**2c**2-a**2D**2-b**2E**2-c**2F**2+DEF),sqrt表示开根号，a**2表示平方
        #D=b**2+c**2-a1**2,E=a**2+c**2-b1**2, F=a**2+b**2-c1**2
        D = b ** 2 + c ** 2 - a1 ** 2
        E = a ** 2 + c ** 2 - b1 ** 2
        F = a ** 2 + b ** 2 - c1 ** 2
        sun=math.sqrt(4*(a ** 2)*(b ** 2)*(c ** 2) - (a ** 2)*(D ** 2) - (b ** 2)*(E ** 2) - (c ** 2)*(F ** 2) + D*E*F)
        V = (1/12)*sun

        """
        # 用线性代数和等体积法求得四面体的体积
        A1 = x2 * (y3 * z4 - y4 * z3) - x3 * (y2 * z4 - y4 * z2) + x4 * (y2 * z3 - y3 * z2)
        A2 = x1 * (y3 * z4 - y4 * z3) - x3 * (y1 * z4 - y4 * z1) + x4 * (y1 * z3 - y3 * z1)
        A3 = x1 * (y2 * z4 - y4 * z2) - x2 * (y1 * z4 - y4 * z1) + x4 * (y1 * z2 - y2 * z1)
        A4 = x1 * (y2 * z3 - y3 * z2) - x2 * (y1 * z3 - y3 * z1) + x3 * (y1 * z2 - y2 * z1)
        v = (A1 - A2 + A3 - A4)/6
        """
        print('四面体体积',round(V,3))
        # 根据等体积法公式1/3*RS=V,求得内切球的半径
        nR = (3*V)/S
        print('内切球半径',round(nR,3))
        narea=4*3.14*nR**2#球体的表面积
        V=(4/3)*3.14*nR**3#球体的体积
        print('内切球表面积',round(narea,3))
        print('内切球体积',round(V,3))
    except:
        print('输入的点不能构成四面体')
def wball(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4):
    #四面体的四个顶点在最小外接球的球面上，因此假设球心为（x,y,z),四个顶点到球心的距离相等，联立四条方程，化简为三元一次方程
    try:
        #化简后，形如：-6.0x-8.0y+2.0z=-26.0', '-4.0x-6.0y-10.0z=-38.0', '-12.0x+0.0y-6.0z=-45.0'
        a=str(2*x1-2*x2)+"x+"+str(2*y1-2*y2)+"y+"+str(2*z1-2*z2)+"z="+str(x1**2-x2**2+y1**2-y2**2+z1**2-z2**2)
        b = str(2 * x1 - 2 * x3) + 'x+'+str(2 * y1 - 2*y3)+ 'y+'+str(2 * z1 - 2 * z3) + 'z=' + str(x1 ** 2 - x3 ** 2 + y1 ** 2 - y3 ** 2 + z1 ** 2 - z3 ** 2)
        c = str(2 * x1 - 2 * x4) + 'x+'+str(2 * y1 - 2 * y4) + 'y+'+str(2 * z1 - 2 * z4) + 'z=' + str(x1 ** 2 - x4 ** 2 + y1 ** 2 - y4 ** 2 + z1 ** 2 - z4 ** 2)
        clist=[a.replace('+-','-'),b.replace('+-','-'),c.replace('+-','-')]
        print('化简后的三元一次方程为：')
        for each in clist:
            print(each)

        # 分别输出三元一次方程的系数矩阵
        a1 = np.array([
            [float(2 * x1 - 2 * x2),float(2*y1-2*y2),float(2*z1-2*z2)],#
            [float(2 * x1 - 2 * x3),float(2 * y1 - 2*y3),float(2 * z1 - 2 * z3)],
            [float(2 * x1 - 2 * x4),float(2 * y1 - 2 * y4),float(2 * z1 - 2 * z4)]
        ])
        # 三元一次方程的值
        b1 = np.array([
            float(x1**2-x2**2+y1**2-y2**2+z1**2-z2**2),
            float(x1 ** 2 - x3 ** 2 + y1 ** 2 - y3 ** 2 + z1 ** 2 - z3 ** 2),
            float(x1 ** 2 - x4 ** 2 + y1 ** 2 - y4 ** 2 + z1 ** 2 - z4 ** 2)])
        # 计算
        x = solve(a1, b1)
        # 打印结果
        print('外接球心坐标为','('+str(round(x[0],3))+','+str(round(x[1],3))+','+str(round(x[2],3))+')')
        R=math.sqrt((x1-x[0])**2+(y1-x[1])**2+(z1-x[2])**2)#根据空间中两点坐标公式求出球体半径
        print('外接球半径',round(R,3))#保留三位小数
        area=4*3.14*R**2#球体的表面积
        print('外接球表面积',round(area,4))
        V=(3/4)*3.14*R**3#球体的体积
        print('外接球体积',round(V,4))
    except:
        print('输入的点不能构成四面体')


def main():
    x1 = float(input('输入第一个点的x值：'))
    y1 = float(input('输入第一个点的y值：'))
    z1 = float(input('输入第一个点的z值：'))

    x2 = float(input('输入第二个点的x值：'))
    y2 = float(input('输入第二个点的y值：'))
    z2 = float(input('输入第二个点的z值：'))

    x3 = float(input('输入第三个点的x值：'))
    y3 = float(input('输入第三个点的y值：'))
    z3 = float(input('输入第三个点的z值：'))

    x4 = float(input('输入第四个点的x值：'))
    y4 = float(input('输入第四个点的y值：'))
    z4 = float(input('输入第四个点的z值：'))
    nball(x1,y1,z1,x2,y2,z2,x3,y3,z3,x4,y4,z4)
    print('+'*14)
    wball(x1, y1, z1, x2, y2, z2, x3, y3, z3, x4, y4, z4)


if __name__=="__main__":
    main()

运行结果
以（6，0，3），（2，3，5），（3，4，-1），（0，0，0）四个点为测试用例.。

结果如下：
输入第一个点的x值：0
输入第一个点的y值：0
输入第一个点的z值：0
输入第二个点的x值：3
输入第二个点的y值：4
输入第二个点的z值：-1
输入第三个点的x值：2
输入第三个点的y值：3
输入第三个点的z值：5
输入第四个点的x值：6
输入第四个点的y值：0
输入第四个点的z值：3
六条边长 [5.0990195135927845, 6.164414002968976, 6.708203932499369, 5.385164807134504, 6.4031242374328485, 6.164414002968976]
四面体表面积 60.609
四面体体积 23.5
内切球半径 1.163
内切球表面积 16.994
内切球体积 6.589
++++++++++++++
化简后的三元一次方程为：
-6.0x-8.0y+2.0z=-26.0
-4.0x-6.0y-10.0z=-38.0
-12.0x+0.0y-6.0z=-45.0
外接球心坐标为 (2.883,1.521,1.734)
外接球半径 3.692
外接球表面积 171.2274
外接球体积 118.5404