每天都能学到自嚎的姿势：
>>什么是牛顿迭代法？
所谓牛顿法(一般称N-R法)是牛顿提出的一种方程求解的方法，过程大概如下：
选取一个初始值(initial guess)x0，过(x0,f(x0))做曲线y=f(x)的切线，找出这条直线与x轴的交点x1，再将x1作为"new guess"，过(x1,f(x1))继续做切线，找出与x轴的交点x2，再将x2作为"new guess"，如此往复，直到f(xn)的值接近我们需要的精度。
用迭代式算法表示：
x_i+1=x_i-f(x_i)/f'(x_i)
>>为什么用牛顿迭代法？
可以证明，如果x0接近方程的根r，则第i次迭代的误差E_i与E_i-1的平方成正比，即二阶收敛。相比之下二分法/黄金分割法/斐波那契法是线性收敛，找到根的速度比牛顿法慢很多。
但这也意味着如果x0远离r，那么需要的迭代次数会很多，甚至无限次迭代也找不到根，这就是一个"poor initial guess"的情况
>>什么是牛顿分形
牛顿分形是Julia集的一种形式，迭代式从z=z^2+c变为z=z-f(z)/f'(z)。考虑一个复平面，从这个平面上任意一点开始，我们会遇到两种情形：
1.从该点经过有限次迭代找到方程的一个根；
2.从该点进行迭代不会收敛，即找不到根。
将这些点在复平面上绘制出来，我们便可以看到一些很奇特的分形。我们得到的分形中，每一个点有两个性质：找到根所需的迭代次数k，以及找到的根r，并根据这些性质进行染色。我今天研究还比较浅，染色原则是k，计划明天加入k-r双重染色的分形。

看图前先进行说明：
>>橙-绿采用的是HSV染色(MATLAB预设的色图)，红色代表最少的迭代次数，绿色代表最多；蓝-青采用JET染色，蓝色最少，青色最多；
>>所有图使用1024x1024矩阵绘制，迭代上限为20次，迭代终止条件为精度<10的5次方。在这种情况下每幅图绘制时间从20-2分钟不等，再精确的对电脑的运算消耗实在太大且画质无明显提升。
首先是最基本的牛顿分形f(x)=x^3-1，图上3个最明显的红点就是方程的解，从r开始到r当然只需要一次迭代，因此次数最少。

以原点为中心放大10倍

将一个“臂”放大25倍可看到无限循环的细节



先简单探索下牛顿分形的性质：
首先考虑方程f(x)=x^3+1，由于解与x^3-1是互补的，不难推测得到的分形与x^3-1对称

再考虑f(x)=x^4-1，由于n次方方程在复数集中有n个解，那么也可推测得到的分形中有n个“峰”，并被分为4半。实际上，在以r为参数染色的分形中，这4半颜色都是不一样的。



探索分形最大的乐趣就在于对迭代方程的变幻可以得到千奇百怪的分形，不怕做不到，只怕想不到。
f(x)=x^8-5x^4-1，虽然形状大变，但红色的“峰”，即解的数量是不变的

f(x)=x^5-3ix^3-2

f(x)=x^(-7+5i)-2

放大8倍

最后是f(x)=sin(x)，为了美观关于直线x=-pi/2放大10倍。我还很脑残地找了e^x，但风扇咆哮了一分多钟后突然想起来，e^x=0不是无解吗-_-b



为了防是不是被不认识的人顶上来骚扰就把源代码放出来吧……转载请事先通知我并在使用的地方注明出处，我不会收取版权费，但我有知情权。
%f为求解的方程，df是导数，使用的时候用funchandler定义
%res是目标分辨率，iter是循环次数，(xc,yc)是图像的中心，xoom是放大倍数
%参数视自己需求增加或减少
function newton(f,df,res,iter,xc,yc,xoom)
%一些乱糟糟的初始化
eol=1e-5;
x0=xc-2.5/xoom;x1=xc+2.5/xoom;
y0=yc-2.5/xoom;y1=yc+2.5/xoom;
x=linspace(x0,x1,res);
y=linspace(y0,y1,res);
[xx,yy]=meshgrid(x,y);
z=xx+yy*1i;
kk=zeros(res,res);
tic
%对每个点进行牛顿迭代
%这个代码并行度很差，因为要对每个点单独进行一次牛顿迭代，执行速度非常慢
%明天的r参数染色我会加入并行化的算法
for m=1:res
for n=1:res
k=0;
t=z(m,n);
ff=z31(t);
while (k<=iter)&&(abs(ff)>eol)
t=t-f(t)/df(t);
ff=f(t);
k=k+1;
end;
kk(m,n)=k;
end;
end;
colormap hsv;
image(x,y,kk);
axis square;
toc
end

完

这次的分形所需的基础知识比较少，有点高中数学知识就能看懂

啊，这次来张高分辨率大倍数的吧，2048x2048 100次迭代，精确到11位，计算用时219秒。

虽然看上去很漂亮 但是完全不知道是什么 = =

ff=z31(t);
这行代码错了

给物理帝跪了！

ff=f(t)
调试的时候用的是z31.m这个函数，最后的.m文件忘改回来了