本文介绍如何用python解方程。
方程定义如下：
输入数据：
年份n
息税前利润y（万元）
总资本x₁=x（万元）
2009
11463.91
279169.98
2010
25899.89
354026.36
2011
20345.11
301854.93
2012
18621.19
297405.49
2013
6178.08
285492.5
2014
7450.71
347817.28
2015
11328.99
389815.43
2016
81077.79
875880.4
2017
24984.51
1113006.97
2018
42230.85
1203954.55
求a, b, c的值
方法1：尝试用excel计算累加数据代入方程，但是发现excel使用浮点数表示数值，存在精度问题。
方法2：我用java的Decimal对象做各种加减乘除，能够保持精度不丢失，但是代码写起来非常费劲，类似(new Decimal().add(new Decimal()))这种。
方法3：stata。但是stata厚厚的文档，大大的学习成本，而我只想帮助mm解一个方程。
方法4：python。都说python好，实际用起来果然名不虚传。
# coding=utf-8
from decimal import *
from sympy import *
if __name__ == '__main__':
# 输入数据x y 单位万元
y = map(Decimal, '11463.91 25899.89 20345.11 18621.19 6178.08 7450.71 11328.99 81077.79 24984.51 42230.85'.split())
# y单位亿元
y = map(lambda n: n / 10000, y)
print 'y = ', y
x = map(Decimal, '279169.98 354026.36 301854.93 297405.49 285492.5 347817.28 389815.43 875880.4 1113006.97 1203954.55'.split())
# x单位亿元
x = map(lambda n: n / 10000, x)
print 'x = ', x
# x1 = x
x1 = x
print 'x1 = ', x1
# x2 = x*x
x2 = map(lambda n: n ** 2, x)
print 'x2 = ', x2
# x1平方
x1_x1 = map(lambda n: n ** 2, x1)
print 'x1_x1 = ', x1_x1
# x2平方
x2_x2 = map(lambda n: n ** 2, x2)
print 'x2_x2 = ', x2_x2
# x1 * x2
x1_x2 = map(lambda n, m: n * m, x1, x2)
print 'x1_x2 = ', x1_x2
# x1 * y
x1_y = map(lambda n, m: n * m, x1, y)
print 'x1_y = ', x1_y
# x2 * y
x2_y = map(lambda n, m: n * m, x2, y)
print 'x2_y = ', x2_y
# sum(y)
sum_y = sum(y)
print 'sum_y = ', sum_y
# sum(x1)
sum_x1 = sum(x1)
print 'sum_x1 ', sum_x1
# sum(x2)
sum_x2 = sum(x2)
print 'sum_x2 ', sum_x2
# sum(x1 * y)
sum_x1_y = sum(x1_y)
print 'sum_x1_y ', sum_x1_y
# sum(x1 * x1)
sum_x1_x1 = sum(x1_x1)
print 'sum_x1_x1 ', sum_x1_x1
# sum(x1 * x2)
sum_x1_x2 = sum(x1_x2)
print 'sum_x1_x2 ', sum_x1_x2
# sum(x2 * y)
sum_x2_y = sum(x2_y)
print 'sum_x2_y ', sum_x2_y
# sum(x2 * x2)
sum_x2_x2 = sum(x2_x2)
print 'sum_x2_x2 ', sum_x2_x2
# 方程1
print sum_y, ' = ', len(y), ' * a + ', sum_x1, ' * b + ', sum_x2, ' * c'
# 方程2
print sum_x1_y, ' = ', sum_x1, ' * a + ', sum_x1_x1, ' * b + ', sum_x1_x2, ' * c'
# 方程3
print sum_x2_y, ' = ', sum_x2, ' * a + ', sum_x1_x2, ' * b + ', sum_x2_x2, ' * c'
# 解方程
a = Symbol('a')
b = Symbol('b')
c = Symbol('c')
print(solve([len(y) * a + sum_x1 * b + sum_x2 * c - sum_y,
sum_x1 * a + sum_x1_x1 * b + sum_x1_x2 * c - sum_x1_y,
sum_x2 * a + sum_x1_x2 * b + sum_x2_x2 * c - sum_x2_y],
[a, b, c]))
计算结果：
y = [Decimal(‘1.146391’), Decimal(‘2.589989’), Decimal(‘2.034511’), Decimal(‘1.862119’), Decimal(‘0.617808’), Decimal(‘0.745071’), Decimal(‘1.132899’), Decimal(‘8.107779’), Decimal(‘2.498451’), Decimal(‘4.223085’)]
x = [Decimal(‘27.916998’), Decimal(‘35.402636’), Decimal(‘30.185493’), Decimal(‘29.740549’), Decimal(‘28.54925’), Decimal(‘34.781728’), Decimal(‘38.981543’), Decimal(‘87.58804’), Decimal(‘111.300697’), Decimal(‘120.395455’)]
x1 = [Decimal(‘27.916998’), Decimal(‘35.402636’), Decimal(‘30.185493’), Decimal(‘29.740549’), Decimal(‘28.54925’), Decimal(‘34.781728’), Decimal(‘38.981543’), Decimal(‘87.58804’), Decimal(‘111.300697’), Decimal(‘120.395455’)]
x2 = [Decimal(‘779.358777332004’), Decimal(‘1253.346635748496’), Decimal(‘911.163987653049’), Decimal(‘884.500254821401’), Decimal(‘815.0596755625’), Decimal(‘1209.768602665984’), Decimal(‘1519.560694660849’), Decimal(‘7671.6647510416’), Decimal(‘12387.845152685809’), Decimal(‘14495.065584657025’)]
x1_x1 = [Decimal(‘779.358777332004’), Decimal(‘1253.346635748496’), Decimal(‘911.163987653049’), Decimal(‘884.500254821401’), Decimal(‘815.0596755625’), Decimal(‘1209.768602665984’), Decimal(‘1519.560694660849’), Decimal(‘7671.6647510416’), Decimal(‘12387.845152685809’), Decimal(‘14495.065584657025’)]
x2_x2 = [Decimal(‘607400.1038044361919084426560’), Decimal(‘1570877.789342073111382166262’), Decimal(‘830219.8123958056305191989964’), Decimal(‘782340.7007791233029464076028’), Decimal(‘664322.27472804775969140625’), Decimal(‘1463540.071996407469712270688’), Decimal(‘2309064.704758161966575133401’), Decimal(‘58854440.05237417450828493056’), Decimal(‘153458707.5269212944961662700’), Decimal(‘210106926.3035085019870868819’)]
x1_x2 = [Decimal(‘21757.357428060001003992’), Decimal(‘44371.774727228591435456’), Decimal(‘27503.934171153197018157’), Decimal(‘26305.523169028362689149’), Decimal(‘23269.342442552703125’), Decimal(‘42077.842480868330340352’), Decimal(‘59234.820560031755710007’), Decimal(‘671946.079080821702464’), Decimal(‘1378775.799822001963708873’), Decimal(‘1745140.016319623543821375’)]
x1_y = [Decimal(‘32.003795254218’), Decimal(‘91.692437811004’), Decimal(‘61.412717548923’), Decimal(‘55.380441363331’), Decimal(‘17.63795504400’), Decimal(‘25.914856862688’), Decimal(‘44.162151083157’), Decimal(‘710.14447136316’), Decimal(‘278.079337720347’), Decimal(‘508.440240078675’)]
x2_y = [Decimal(‘893.449888104413397564’), Decimal(‘3246.153999775611406544’), Decimal(‘1853.773155683992374039’), Decimal(‘1647.044730007772408719’), Decimal(‘503.5503880399170000’), Decimal(‘901.363502556947364864’), Decimal(‘1721.508791420581171251’), Decimal(‘62200.1623635353126064’), Decimal(‘30950.424109573012181859’), Decimal(‘61213.894044581312422125’)]
sum_y = 24.958103
sum_x1 544.842389
sum_x2 41927.334116828717
sum_x1_y 1824.868404129503
sum_x1_x1 41927.334116828717
sum_x1_x2 4040382.490201370151316361
sum_x2_y 165131.324973278872333365
sum_x2_x2 430647839.3406080264242731083
24.958103 = 10 * a + 544.842389 * b + 41927.334116828717 * c
1824.868404129503 = 544.842389 * a + 41927.334116828717 * b + 4040382.490201370151316361 * c
165131.324973278872333365 = 41927.334116828717 * a + 4040382.490201370151316361 * b + 430647839.3406080264242731083 * c
{c: -0.00208988296109613, b: 0.337762078993800, a: -7.14457738576440}