原本是想直接续在这帖:http://tieba.baidu.com/p/2134509239的后面直接回复的,但是考虑到那帖楼层数已经较多,默默发在那边无法满足我的虚荣心(殴)所以决定重开一帖。(其实@隨意超 似乎已经在那帖11楼找对了方向,但却似乎没有引起后续的参与者的重视。)
先上代码和图再讲解吧。2GB内存2GHz老爷机,3 min出图:
modifiedroots[c_List] :=
Block[{a = DiagonalMatrix[ConstantArray[1., Length@c - 1], -1]},
a[[1]] = -c;
Eigenvalues[a]];
s = 1000;
l = ConstantArray[0., {s + 1, s + 1}];
l[[#, #2]] += 1. & @@@
Round[1 + s Rescale@
Function[z, {Im@z, Re@z}, Listable]@
Flatten[modifiedroots /@ Tuples[{-1., 1.}, 18]]]; // AbsoluteTiming
ArrayPlot[UnitStep[80 - l] l, ColorFunction -> "AvocadoColors"]
当然,我说“复现”其实是标题党了,因为原图算到了24次多项式,我这幅图只算到了19次。(这是2G内存的极限,次数每加1,内存占用大约要翻一番。)但是如大家所见,这幅图的效果已经和那幅Sam Derbyshire画了3天3夜的图相差无几了。
下面简单讲解一下这段代码。要做到在多数人可以接受的时间内使用多数人可以接受的内存获得和那幅图相似的效果,需要:
1 省略对这个特别的问题而言其实不必要的(主要是符号)计算。
2 昨晚对像素点十分有限的屏幕而言其实不必要的仅有着细微差别的数据。
3 正确地注意到原图的颜色代表的是什么。
先说第一点。因为这幅图表示的是幂方程的根,所以我们自然需要求解幂方程。如何解方程?我们的第一反应是Solve。是的,Solve可以求解这个问题,但是,作为一个通用求解器,它对这个问题而言太慢了。有没有更为专用的函数可以求解这个问题?在本文开头的帖子里@草红样 找到了NRoots,一个专解多项式方程的函数,从而大大提速了求解。可是,这是否是终点?不,NRoots依旧需要进行多项式结构的分析,我们能不能把这一部也给省掉,直接使用(以确定次序排列的)多项式系数表作为输入来求根?再一次,答案是肯定的,那就是我们上面使用的modifiedroots。(对于NRoots的提速问题的讨论其实尚未结束,有兴趣的参考这帖:http://tieba.baidu.com/p/3575212088)
另,不难注意到“全体系数为正负1的全体n次幂方程”其实有一半是重复的。上面的modifiedroots已经利用了这个特性。
再说后两点。(这两点其实是相联系的,所以一起讲。)对方程的求解将会产生大量的数据,这倒是其次,真正够呛的是如何将这些数据绘成图形。(具体参看本文开头帖子中大伙儿的哀嚎。)正常情况下,要画这么多的点很花时间,还吃内存,更麻烦的是,还怎么画也不像Sam Derbyshire那幅“震撼”的图!要跨过这个槛,需要另一个技巧,好吧,其实这个技巧我是从这帖学来的:http://mathematica.stackexchange.com/questions/50839/smooth-peter-de-jong-attractor#comment149926_50839 (注意下面的答案最好也看看),这里的重点是,Sam那幅图中的颜色明暗到底表现的是什么?复平面上同一点处重合的点量?我们不难想像完全重合的点其实是很少的,是的,颜色的明暗其实表现的只是某一个点附近小区域内的点量。如何决定这个小区域?对于这点,我们的计算机屏幕其实已经帮我们限制好了:屏幕上的像素点是很有限的。于是这又反过来启发了我们:既然最终在屏幕上呈现的点数有限,我们又何必把所有的点都拿去画图?只要事先统计好每个像素内有多少点,再把这个统计结果拿去画图就好了嘛。哪个函数可以经济地画出这种原始数据点数极多的明暗图?我们有ArrayPlot。ArrayPlot可以调配色方案嗯。
另,最后画图时使用UnitStep剪掉了实轴上的点,因为影响视觉效果。
感觉其实还有些地方地方可以解说下不过现在也累了所以算了。
最后,其实即便抛开方程阶数,单论配色方案,我也没有完全“复现”人家的图。(小声:不过我自以为我选的方案更好看。)有谁有兴趣可以研究下怎么用Blend把原图的颜色混合出来。(应该不是某个Mathematica里现成的配色方案吧……)
先上代码和图再讲解吧。2GB内存2GHz老爷机,3 min出图:
modifiedroots[c_List] :=
Block[{a = DiagonalMatrix[ConstantArray[1., Length@c - 1], -1]},
a[[1]] = -c;
Eigenvalues[a]];
s = 1000;
l = ConstantArray[0., {s + 1, s + 1}];
l[[#, #2]] += 1. & @@@
Round[1 + s Rescale@
Function[z, {Im@z, Re@z}, Listable]@
Flatten[modifiedroots /@ Tuples[{-1., 1.}, 18]]]; // AbsoluteTiming
ArrayPlot[UnitStep[80 - l] l, ColorFunction -> "AvocadoColors"]
当然,我说“复现”其实是标题党了,因为原图算到了24次多项式,我这幅图只算到了19次。(这是2G内存的极限,次数每加1,内存占用大约要翻一番。)但是如大家所见,这幅图的效果已经和那幅Sam Derbyshire画了3天3夜的图相差无几了。
下面简单讲解一下这段代码。要做到在多数人可以接受的时间内使用多数人可以接受的内存获得和那幅图相似的效果,需要:
1 省略对这个特别的问题而言其实不必要的(主要是符号)计算。
2 昨晚对像素点十分有限的屏幕而言其实不必要的仅有着细微差别的数据。
3 正确地注意到原图的颜色代表的是什么。
先说第一点。因为这幅图表示的是幂方程的根,所以我们自然需要求解幂方程。如何解方程?我们的第一反应是Solve。是的,Solve可以求解这个问题,但是,作为一个通用求解器,它对这个问题而言太慢了。有没有更为专用的函数可以求解这个问题?在本文开头的帖子里@草红样 找到了NRoots,一个专解多项式方程的函数,从而大大提速了求解。可是,这是否是终点?不,NRoots依旧需要进行多项式结构的分析,我们能不能把这一部也给省掉,直接使用(以确定次序排列的)多项式系数表作为输入来求根?再一次,答案是肯定的,那就是我们上面使用的modifiedroots。(对于NRoots的提速问题的讨论其实尚未结束,有兴趣的参考这帖:http://tieba.baidu.com/p/3575212088)
另,不难注意到“全体系数为正负1的全体n次幂方程”其实有一半是重复的。上面的modifiedroots已经利用了这个特性。
再说后两点。(这两点其实是相联系的,所以一起讲。)对方程的求解将会产生大量的数据,这倒是其次,真正够呛的是如何将这些数据绘成图形。(具体参看本文开头帖子中大伙儿的哀嚎。)正常情况下,要画这么多的点很花时间,还吃内存,更麻烦的是,还怎么画也不像Sam Derbyshire那幅“震撼”的图!要跨过这个槛,需要另一个技巧,好吧,其实这个技巧我是从这帖学来的:http://mathematica.stackexchange.com/questions/50839/smooth-peter-de-jong-attractor#comment149926_50839 (注意下面的答案最好也看看),这里的重点是,Sam那幅图中的颜色明暗到底表现的是什么?复平面上同一点处重合的点量?我们不难想像完全重合的点其实是很少的,是的,颜色的明暗其实表现的只是某一个点附近小区域内的点量。如何决定这个小区域?对于这点,我们的计算机屏幕其实已经帮我们限制好了:屏幕上的像素点是很有限的。于是这又反过来启发了我们:既然最终在屏幕上呈现的点数有限,我们又何必把所有的点都拿去画图?只要事先统计好每个像素内有多少点,再把这个统计结果拿去画图就好了嘛。哪个函数可以经济地画出这种原始数据点数极多的明暗图?我们有ArrayPlot。ArrayPlot可以调配色方案嗯。
另,最后画图时使用UnitStep剪掉了实轴上的点,因为影响视觉效果。
感觉其实还有些地方地方可以解说下不过现在也累了所以算了。
最后,其实即便抛开方程阶数,单论配色方案,我也没有完全“复现”人家的图。(小声:不过我自以为我选的方案更好看。)有谁有兴趣可以研究下怎么用Blend把原图的颜色混合出来。(应该不是某个Mathematica里现成的配色方案吧……)