数学应用软件大型实验实验报告：涂色问题

数学应用软件大型实验实验报告 实验序号： 日期：2012 年6 月 21 日 班级 姓名 学号 实验 名称
数学软件实验：通过牛顿迭代公式得出彩色图形 问题背景描述：
牛顿迭代公式可以直接用来求解复数方程Z^3-1=0，在复平面上他的三个根 是Z(1)=1，Z(2）＝-1/2+1/2*i*3^(1/2)，Z(3）＝-1/2－1/2*i*3^(1/2) 选择中心位于坐标原点边长为2的正方形内的点为初始值，把收敛到三个不同
根的初始值分别标上不同的颜色。

只要计算足够多的点，你将得到牛顿法收敛域的彩色图形。

实验目的：
1.了解牛顿迭代法解复数方程时收敛域的结构。

2.通过实验巩固对使用MATLAB 绘制图形等知识加深学习。

实验原理与数学模型：
牛顿迭代法可以用于求解复数方程01=-n z ，例如对016
=-z ，该方程在复平面上有六个根，选择中心位于坐标原点，边长为4的正方形内的任意点作初始值进行迭代，把收敛到六个根的初值归为另外六类，分别以不同颜色做图。

对充分多的初始点进行实验，绘出牛顿迭代法对该方程的收敛域彩色图。

绘伪彩色图命令pcolor()
该命令主要用于绘制矩阵色图，根据矩阵中元素数据的大小不同绘不同颜色。

常常与命令shading interp 结合使用效果会更好。

求多项式零点命令roots()
该命令用于求多项式P(x) = a 1x n + a 2 x (n-1) + ... + a n x + a n+1 的全部零点。

实验所用软件及版本：
MATLAB 7.0.1
主要内容（要点）：
记1)(3-=z z f ，则 23)(z z f ='，所以牛顿迭代公式为3/)/1(21n n n n z z z z --=+
取迭代次数为10。

考虑以坐标原点为中心的正方形区域2}y 2- 2, x 2- | y) {(x,≤≤≤≤=Ω 取步长h=0.02，在区域内取离散网格点
)200,...,1,0,(,22=⎩⎨⎧⨯+-=⨯+-=k j h k y h j x k
j 由此可以构造出规则的复数)200,...,1,0,(,=+=k j iy x z k j jk 。

对于这些点，逐一用牛顿迭代法取初值进行迭代实验。

为了记录实验结果，构造三个阶数均为201×201 矩阵：Z 1、Z 2、Z 3，开始时这三个矩阵都设为全零矩阵。

如果以jk z 为初值的迭代实验结果是收敛到第一个根，则将Z 1的第j 行第k 列的元素改写为1；如果以jk z 为初值的迭代实验结果是收敛到第二个根，则将Z 2的第j 行第k 列的元素改写为1；如果以jk z 为初值的迭代实验结果是收敛到第三个根，则将Z 3的第j 行第k 列的元素改写为1。

令Z = Z 1 + 2Z 2 + 3Z 3，则矩阵Z 的元素由“1”、“2”、“3”这三个元素组成。

位置上数据为“1”，说明这一复数做初值导致牛顿迭代法收敛到第一个根；位置上数据为“2”，说明这一复数收敛到第二个根；位置上数据为“3”，说明这一复数收敛到第三个根。

将这一矩阵用MATLAB 作图命令pcolor()作用，将绘出收敛区域色图。

实验过程记录（含：基本步骤、主要程序清单及异常情况记录等）：
程序过程：
X=roots([1,0,0,-1]);
r1=X(1);r2=X(2);r3=X(3);
h=0.02;N=1+4/h;
z0(N,N)=0;z1=z0;z2=z0;z3=z0;
t=(-2:h:2)+eps;
[x,y]=meshgrid(t);
z=x+y*i;
for j=1:N
for k=1:N
p=z(j,k);
for n=1:10
p=p-(p-1/p^2)/3;
end
if abs(p-r1)<0.01
z1(j,k)=1;
elseif abs(p-r2)<0.01
z2(j,k)=1;
elseif abs(p-r3)<0.01
z3(j,k)=1;
end
end
end
Z=z1+2*z2+3*z3;
figure(1)
pcolor(x,y,Z),shading interp
实验结果报告与实验总结：
通过以上过程得出结果如图（figure 1)：
2
1.5
1
0.5
-0.5
-1
-1.5
-1.5-1-0.500.51 1.52
-1.5-1-0.500.51 1.52-1.5
-1-0.5
00.511.5
2
思考与深入：
程序中已经涉及到了另一部分，即当不收敛时的彩色图形。

加入不收敛的部分时，程序如下：
X=roots([1,0,0,-1]); r1=X(1);r2=X(2);r3=X(3);h=0.02;N=1+4/h; z0(N,N)=0;z1=z0;z2=z0;z3=z0;
t=(-2:h:2)+eps;[x,y]=meshgrid(t); z=x+y*i;for j=1:N
for k=1:N
p=z(j,k); for n=1:10
p=p-(p-1/p^2)/3;
end
if abs(p-r1)<0.01
z1(j,k)=1;
elseif abs(p-r2)<0.01
z2(j,k)=1; elseif abs(p-r3)<0.01 z3(j,k)=1; else
z0(j,k)=1;
end
end end Z=z0+2*z1+3*z2+4*z3;figure(1)
pcolor(x,y,Z),shading interp figure(2) pcolor(x,y,z0),shading interp
上图为收敛时（figure 1),下图为不收敛时（figure 2)。

-1.5-1-0.500.51 1.52
-1.5
-1
-0.5
0.5
1
1.5
2。