不动点迭代法和牛顿法非线性方程组求解.docx

不动点迭代法和牛顿法非线性方程组求解.docx

《不动点迭代法和牛顿法非线性方程组求解.docx》由会员分享，可在线阅读，更多相关《不动点迭代法和牛顿法非线性方程组求解.docx（17页珍藏版）》请在冰豆网上搜索。

不动点迭代法和牛顿法非线性方程组求解

《MATLAB程序设计实践》课程考核

１、编程实现以下科学计算算法，并举一例应用之。

（参考书籍《精

通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９

年）

“不动点迭代法和牛顿法非线性方程组求解”

（1）.不动点迭代法非线性方程组求解

（a）.算法说明：

设含有n个未知数与n个方程的非线性方程组记为：

F（x）=0，然后把上述方程组改为便于迭代的等价形式：

x=φ（x），由此就可以构造不动点迭代法的迭代公式：

如果得到的序列{xk}满足

，则

就是φ的不动点，这样就可以求出非线性方程组的解。

在MATLAB中编程实现的非线性方程组的不动点迭代法的函数为：

mulStablePoint。

功能：

用不动点迭代法求非线性方程组的一个解。

调用格式：

[r,n]=mulStablePoint（x0,eps）。

其中，x0为初始迭代向量；

eps为迭代精度；

r为求出的解向量；

n为迭代步数。

（b）.流程图：

（c）.源程序代码：

function[r,n]=mulStablePoint（x0,eps）

%不动点迭代法求非线性方程组的一个解

%初始迭代向量：

x0

%迭代精度：

eps

%解向量：

r

%迭代步数：

n

ifnargin==1

eps=1.0e-4;

end

r=myf（x0）;

n=1;

tol=1;

whiletol>eps

x0=r;

r=myf（x0）;%迭代公式

tol=norm（r-x0）;%注意矩阵的误差求法，norm为矩阵的欧几里德范数

n=n+1;

if（n>100000）%迭代步数控制

disp（'迭代步数太多，可能不收敛！

'）;

return;

end

end

举例说明：

解：

首先建立myf.m函数文件，输入以下内容：

functionf=myf（x）

f

（1）=0.5*sin（x

（1））+0.1*cos（x

（2）*x

（1））-x

（1）;

f

（2）=0.5*cos（x

（1））-0.1*sin（x

（2））-x

（2）;

在MATLAB命令窗口中输入：

（2）.牛顿法非线性方程组求解

（a）.算法说明：

设非线性方程组为

，其中

，

牛顿迭代法的迭代公式为：

求解步骤为：

（1）给出初始值

；

（2）对n=1,2,3…计算F（xn）和F’（xn）；

（3）求出xn+1，并进行精度控制。

更一般的牛顿法迭代公式为：

当

=F’（x0）时，就得到简化牛顿法。

在MATLAB中编程实现的非线性方程组的牛顿迭代法的函数为：

mulNewton。

功能：

用牛顿迭代法求非线性方程组的一个解。

调用格式：

[r,n]=mulNewton（x0,eps）。

其中，x0为初始迭代向量；

eps为迭代精度；

r为求出的解向量；

n为迭代步数。

（b）.流程图：

（c）.源程序代码：

function[r,n]=mulNewton（x0,eps）

%牛顿迭代法求非线性方程组的一个解

%初始迭代向量x0

%迭代精度eps

%解向量r

%迭代步数n

ifnargin==1

eps=1.0e-4;

end

r=x0-myf（x0）/dmyf（x0）;

n=1;

tol=1;

whiletol>eps

x0=r;

r=x0-myf（x0）/dmyf（x0）;%核心迭代公式

tol=norm（r-x0）;

n=n+1;

if（n>100000）%迭代步数控制

disp（'迭代步数太多，可能不收敛！

'）;

return;

end

end

另一种方法为简化牛顿迭代法，如下：

在MATLAB中编程实现的非线性方程组的简化牛顿迭代法的函数为：

mulSimNewton。

功能：

用简化牛顿迭代法求非线性方程组的一个解。

调用格式：

[r,n]=mulSimNewton（x0,eps）。

其中，x0为初始迭代向量；

eps为迭代精度；

r为求出的解向量；

n为迭代步数。

源程序代码：

function[r,n]=mulSimNewton（x0,eps）

%简化牛顿迭代法求非线性方程组的一个解

%初始迭代向量x0

%迭代精度eps

%解向量r

%迭代步数n

ifnargin==1

eps=1.0e-4;

end

r=x0-myf（x0）/dmyf（x0）;

c=dmyf（x0）;

n=1;

tol=1;

whiletol>eps

x0=r;

r=x0-myf（x0）/c;

tol=norm（r-x0）;

n=n+1;

if（n>100000）

disp（'‘迭代步数太多，可能不收敛！

'）;

return;

end

end

举例说明：

解：

首先建立myf.m函数文件，输入以下内容：

functionf=myf（x）

f

（1）=0.5*sin（x

（1））+0.1*cos（x

（2）*x

（1））-x

（1）;

f

（2）=0.5*cos（x

（1））-0.1*sin（x

（2））-x

（2）;

f=[f

（1）f

（2）];

再建立dmyf.m导数的雅可比矩阵，输入以下内容：

functiondf=dmyf（x）

df=[0.5*cos（x

（1））-0.1*x

（2）*sin（x

（2）*x

（1））-1-0.1*x

（1）*sin（x

（2）*x

（1））-0.5*sin（x

（1））-0.1*cos（x

（2））-1];

然后在MATLAB命令窗口中输入：

2、编程解决以下科学计算问题

1）

5．有3个多项式

试进行下列操作：

（1）求

＝

＋

（2）求

的根

（3）当x取矩阵A的每一个元素，求

的值。

其中：

A＝

（4）当以矩阵A为自变量时，求

的值。

其中A的值与第（3）题相同。

流程图

源程序代码：

%求积函数conv（x,y）

%求平方根函数roots（p）

%逐一取用矩阵中的数值函数polyval（p,x）

%取用矩阵的函数polyvalm（p,A）是按照矩阵运算规则计算多项式的值

%多项式p2，p3

%矩阵A

p1=[12405];

p2=[12];

p3=[123];

p0=conv（p2,p3）%对p2和p3求积

p=p1+[0p0]%对p1和p0进行求和

x=roots（p）%对p进行求根

x=[-11.2-1.4;0.7523.5;052.5];

x0=polyval（p,x）%将矩阵x中的每一个数值代入p中

A=[-11.2-1.4;0.7523.5;052.5];

x0=polyvalm（p,A）%将矩阵A代入p中

结果

2）

2.用三次多项式拟合下面数据，做出图形。

x=[00.20.40.60.81]

y=[07.7810.688.373.970]

解：

【.m文件程序代码】

>>x=0:

0.2:

1;

>>y=[07.7810.688.373.970];

>>a=polyfit（x,y,3）

a=

41.5625-101.607160.0768-0.1179

>>x1=[0:

0.05:

1];

>>y1=a（4）+a（3）*x1+a

（2）*x1.^2+a

（1）*x1.^3;

>>plot（x,y,'*'）

>>holdon

>>plot（x1,y1,'-r'）

>>

运行结果：

流程图：

3.拟合函数有如下形式：

y=αexp（βx）

试确定系数，并分别用线性尺度和对数尺度做出拟合曲线的图形。

x=[0.01290.02470.05300.15500.30100.47100.80201.27001.43002.4600]

y=[9.56008.18455.26122.79172.26111.73401.23701.06741.11710.7620]

解：

【.m文件程序代码】

functionjisuan

x=[0.01290.02470.0530.1550.3010.4710.8021.271.432.46];

y=[9.568.18455.26162.79172.26111.7341.2371.06741.11710.762];

p=polyfit（x,log（y）,1）;

a=exp（p

（2））

b=p

（1）

x0=0.0129:

0.0001:

2.46;

y0=4.4717.*exp（-0.9238.*x0）;

subplot（2,1,1）

plot（x0,y0）

holdon

subplot（2,1,2）

loglog（x0,y0）

holdon

>>jisuan

a=

4.4717

b=

-0.9238

流程图：