Python 求超越方程解

2016-05-26 19:20:20 +08:00
 asxaqz
现在我有若干带参数的超越方程,比如:

我希望把每一个超越方程的解表示成参数的函数的形式(不是显式解,因为方程中中有一个比较长的幂级数, sympy 我觉得是做不了的)。我现在的想法是对参数空间进行搜索,把参数空间的每一个点对应的解求出来,然后进行插值。这些用 scipy 都有函数可以做。

但是我担心最后精度和稳定性的问题,因为我后面的超越方程要用到前面超越方程的解。不知道有没有人做过相关的事情?或者有现成的包推荐吗?谢谢!
6076 次点击
所在节点    Python
8 条回复
ayaseangle
2016-05-26 19:29:55 +08:00
为什么规定 python
xiaket
2016-05-26 19:50:14 +08:00
话说, 写 python 去调 Mathematica 网站的 API 算不算违规...
strider0505
2016-05-26 20:08:05 +08:00
#!/usr/bin/python3
from math import log, exp
## ln(f)=ln(p)+q1*f=ln(p)+q1*exp[ln(f)]
def f(p,q):
lp=log(p)
epsilon=1e-5
## initial guess
y0=0
y1=1
while (y1-y0)>epsilon or (y0-y1)>epsilon:
y1=lp+q*exp(y1)
y0=y1
return exp(y1)

if __name__=="__main__":
for p in range(10):
for q in range(10):
print("p=%s, q=%s, f=%s", p/10+0.1, q/10+0.1, f(p/10+0.1,q/10+0.1))



结果:
p=%s, q=%s, f=%s 0.1 0.1 0.13123614957214594
p=%s, q=%s, f=%s 0.1 0.2 0.17222926954522666
p=%s, q=%s, f=%s 0.1 0.30000000000000004 0.22602706178738802
p=%s, q=%s, f=%s 0.1 0.4 0.2966292128808233
p=%s, q=%s, f=%s 0.1 0.5 0.38928475749095637
p=%s, q=%s, f=%s 0.1 0.6 0.510882326602397
p=%s, q=%s, f=%s 0.1 0.7 0.6704622942775809
p=%s, q=%s, f=%s 0.1 0.7999999999999999 0.8798888993429671
p=%s, q=%s, f=%s 0.1 0.9 1.1547323120104451
p=%s, q=%s, f=%s 0.1 1.0 1.5154262241479266
p=%s, q=%s, f=%s 0.2 0.1 0.26247229914429193
p=%s, q=%s, f=%s 0.2 0.2 0.3444585390904532
p=%s, q=%s, f=%s 0.2 0.30000000000000004 0.45205412357477603
p=%s, q=%s, f=%s 0.2 0.4 0.5932584257616464
p=%s, q=%s, f=%s 0.2 0.5 0.7785695149819126
p=%s, q=%s, f=%s 0.2 0.6 1.0217646532047937
p=%s, q=%s, f=%s 0.2 0.7 1.3409245885551615
p=%s, q=%s, f=%s 0.2 0.7999999999999999 1.759777798685934
p=%s, q=%s, f=%s 0.2 0.9 2.3094646240208903
p=%s, q=%s, f=%s 0.2 1.0 3.0308524482958528
p=%s, q=%s, f=%s 0.30000000000000004 0.1 0.3937084487164379
p=%s, q=%s, f=%s 0.30000000000000004 0.2 0.5166878086356799
p=%s, q=%s, f=%s 0.30000000000000004 0.30000000000000004 0.678081185362164
p=%s, q=%s, f=%s 0.30000000000000004 0.4 0.8898876386424697
p=%s, q=%s, f=%s 0.30000000000000004 0.5 1.1678542724728689
p=%s, q=%s, f=%s 0.30000000000000004 0.6 1.5326469798071907
p=%s, q=%s, f=%s 0.30000000000000004 0.7 2.0113868828327424
p=%s, q=%s, f=%s 0.30000000000000004 0.7999999999999999 2.639666698028901
p=%s, q=%s, f=%s 0.30000000000000004 0.9 3.4641969360313354
p=%s, q=%s, f=%s 0.30000000000000004 1.0 4.546278672443779
p=%s, q=%s, f=%s 0.4 0.1 0.5249445982885838
p=%s, q=%s, f=%s 0.4 0.2 0.6889170781809064
p=%s, q=%s, f=%s 0.4 0.30000000000000004 0.9041082471495521
p=%s, q=%s, f=%s 0.4 0.4 1.186516851523293
p=%s, q=%s, f=%s 0.4 0.5 1.557139029963825
p=%s, q=%s, f=%s 0.4 0.6 2.0435293064095874
p=%s, q=%s, f=%s 0.4 0.7 2.681849177110323
p=%s, q=%s, f=%s 0.4 0.7999999999999999 3.5195555973718675
p=%s, q=%s, f=%s 0.4 0.9 4.61892924804178
p=%s, q=%s, f=%s 0.4 1.0 6.061704896591705
p=%s, q=%s, f=%s 0.5 0.1 0.6561807478607297
p=%s, q=%s, f=%s 0.5 0.2 0.8611463477261331
p=%s, q=%s, f=%s 0.5 0.30000000000000004 1.13013530893694
p=%s, q=%s, f=%s 0.5 0.4 1.4831460644041161
p=%s, q=%s, f=%s 0.5 0.5 1.9464237874547814
p=%s, q=%s, f=%s 0.5 0.6 2.554411633011984
p=%s, q=%s, f=%s 0.5 0.7 3.352311471387903
p=%s, q=%s, f=%s 0.5 0.7999999999999999 4.399444496714834
p=%s, q=%s, f=%s 0.5 0.9 5.773661560052224
p=%s, q=%s, f=%s 0.5 1.0 7.57713112073963
p=%s, q=%s, f=%s 0.6 0.1 0.7874168974328756
p=%s, q=%s, f=%s 0.6 0.2 1.0333756172713595
p=%s, q=%s, f=%s 0.6 0.30000000000000004 1.3561623707243278
p=%s, q=%s, f=%s 0.6 0.4 1.7797752772849391
p=%s, q=%s, f=%s 0.6 0.5 2.3357085449457373
p=%s, q=%s, f=%s 0.6 0.6 3.065293959614381
p=%s, q=%s, f=%s 0.6 0.7 4.022773765665484
p=%s, q=%s, f=%s 0.6 0.7999999999999999 5.279333396057801
p=%s, q=%s, f=%s 0.6 0.9 6.92839387206267
p=%s, q=%s, f=%s 0.6 1.0 9.092557344887554
p=%s, q=%s, f=%s 0.7 0.1 0.9186530470050215
p=%s, q=%s, f=%s 0.7 0.2 1.205604886816586
p=%s, q=%s, f=%s 0.7 0.30000000000000004 1.5821894325117158
p=%s, q=%s, f=%s 0.7 0.4 2.0764044901657623
p=%s, q=%s, f=%s 0.7 0.5 2.7249933024366935
p=%s, q=%s, f=%s 0.7 0.6 3.5761762862167776
p=%s, q=%s, f=%s 0.7 0.7 4.6932360599430645
p=%s, q=%s, f=%s 0.7 0.7999999999999999 6.159222295400768
p=%s, q=%s, f=%s 0.7 0.9 8.083126184073114
p=%s, q=%s, f=%s 0.7 1.0 10.607983569035483
p=%s, q=%s, f=%s 0.7999999999999999 0.1 1.0498891965771675
p=%s, q=%s, f=%s 0.7999999999999999 0.2 1.3778341563618128
p=%s, q=%s, f=%s 0.7999999999999999 0.30000000000000004 1.8082164942991037
p=%s, q=%s, f=%s 0.7999999999999999 0.4 2.3730337030465853
p=%s, q=%s, f=%s 0.7999999999999999 0.5 3.1142780599276496
p=%s, q=%s, f=%s 0.7999999999999999 0.6 4.087058612819174
p=%s, q=%s, f=%s 0.7999999999999999 0.7 5.363698354220644
p=%s, q=%s, f=%s 0.7999999999999999 0.7999999999999999 7.039111194743734
p=%s, q=%s, f=%s 0.7999999999999999 0.9 9.23785849608356
p=%s, q=%s, f=%s 0.7999999999999999 1.0 12.123409793183411
p=%s, q=%s, f=%s 0.9 0.1 1.1811253461493134
p=%s, q=%s, f=%s 0.9 0.2 1.5500634259070394
p=%s, q=%s, f=%s 0.9 0.30000000000000004 2.034243556086492
p=%s, q=%s, f=%s 0.9 0.4 2.669662915927409
p=%s, q=%s, f=%s 0.9 0.5 3.503562817418606
p=%s, q=%s, f=%s 0.9 0.6 4.597940939421571
p=%s, q=%s, f=%s 0.9 0.7 6.034160648498226
p=%s, q=%s, f=%s 0.9 0.7999999999999999 7.919000094086702
p=%s, q=%s, f=%s 0.9 0.9 10.392590808094004
p=%s, q=%s, f=%s 0.9 1.0 13.638836017331336
p=%s, q=%s, f=%s 1.0 0.1 1.3123614957214593
p=%s, q=%s, f=%s 1.0 0.2 1.722292695452266
p=%s, q=%s, f=%s 1.0 0.30000000000000004 2.26027061787388
p=%s, q=%s, f=%s 1.0 0.4 2.966292128808232
p=%s, q=%s, f=%s 1.0 0.5 3.8928475749095623
p=%s, q=%s, f=%s 1.0 0.6 5.108823266023968
p=%s, q=%s, f=%s 1.0 0.7 6.704622942775806
p=%s, q=%s, f=%s 1.0 0.7999999999999999 8.798888993429669
p=%s, q=%s, f=%s 1.0 0.9 11.54732312010445
p=%s, q=%s, f=%s 1.0 1.0 15.154262241479262
w466397352
2016-05-26 20:18:07 +08:00
数值分析啊…可惜当初全还给老师了_(:з」∠)_
debiann
2016-05-26 20:20:42 +08:00
用 mathematica 吧
Mutoo
2016-05-26 20:51:15 +08:00
@strider0505
>>> "value=%s"%(1.0)
'value=1.0'
thekoc
2016-05-26 22:17:07 +08:00
不知道 sympy 行不行
thekoc
2016-05-26 22:17:51 +08:00
上面那条画一个删除线,还是 mathematica 吧= =

这是一个专为移动设备优化的页面(即为了让你能够在 Google 搜索结果里秒开这个页面),如果你希望参与 V2EX 社区的讨论,你可以继续到 V2EX 上打开本讨论主题的完整版本。

https://www.v2ex.com/t/281504

V2EX 是创意工作者们的社区,是一个分享自己正在做的有趣事物、交流想法,可以遇见新朋友甚至新机会的地方。

V2EX is a community of developers, designers and creative people.

© 2021 V2EX