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

《MATLAB程序设计实践》课程考核１、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“不动点迭代法和牛顿法非线性方程组求解”(1).不动点迭代法非线性方程组求解(a).算法说明：设含有n个未知数与n个方程的非线性方程组记为：F(x)=0，然后把上述方程组改为便于迭代的等价形式：x=φ(x)，由此就可以构造不动点迭代法的迭代公式：如果得到的序列{xk}满足错误!未找到引用源。，则错误!未找到引用源。就是φ的不动点，这样就可以求出非线性方程组的解。在MATLAB中编程实现的非线性方程组的不动点迭代法的函数为：mulStablePoint。功能：用不动点迭代法求非线性方程组的一个解。调用格式：[r,n]=mulStablePoint(x0,eps)。其中，x0为初始迭代向量；eps为迭代精度；r为求出的解向量；n为迭代步数。(b).流程图：(c).源程序代码：NYYNNY开始输入初始迭代值nargin==1r=myf(x0)n=1tol=1tolepsx0=rr=myf(x0)tol=norm(r-x0)n=n+1n100000输出结果eps=1.0e-4迭代步数太多，可能不收敛！结束function[r,n]=mulStablePoint(x0,eps)%不动点迭代法求非线性方程组的一个解%初始迭代向量：x0%迭代精度：eps%解向量：r%迭代步数：nifnargin==1eps=1.0e-4;endr=myf(x0);n=1;tol=1;whiletolepsx0=r;r=myf(x0);%迭代公式tol=norm(r-x0);%注意矩阵的误差求法，norm为矩阵的欧几里德范数n=n+1;if(n100000)%迭代步数控制disp('迭代步数太多，可能不收敛！');return;endend举例说明：解：首先建立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)%牛顿迭代法求非线性方程组的一个解NYYNNY开始输入初始迭代值nargin==1r=x0-myf(x0)/dmyf(x0)n=1tol=1tolepsx0=rr=x0-myf(x0)/dmyf(x0)tol=norm(r-x0)n=n+1n100000输出结果eps=1.0e-4迭代步数太多，可能不收敛！结束%初始迭代向量x0%迭代精度eps%解向量r%迭代步数nifnargin==1eps=1.0e-4;endr=x0-myf(x0)/dmyf(x0);n=1;tol=1;whiletolepsx0=r;r=x0-myf(x0)/dmyf(x0);%核心迭代公式tol=norm(r-x0);n=n+1;if(n100000)%迭代步数控制disp('迭代步数太多，可能不收敛！');return;endend另一种方法为简化牛顿迭代法，如下：在MATLAB中编程实现的非线性方程组的简化牛顿迭代法的函数为：mulSimNewton。功能：用简化牛顿迭代法求非线性方程组的一个解。调用格式：[r,n]=mulSimNewton(x0,eps)。其中，x0为初始迭代向量；eps为迭代精度；r为求出的解向量；n为迭代步数。源程序代码：function[r,n]=mulSimNewton(x0,eps)%简化牛顿迭代法求非线性方程组的一个解%初始迭代向量x0%迭代精度eps%解向量r%迭代步数nifnargin==1eps=1.0e-4;endr=x0-myf(x0)/dmyf(x0);c=dmyf(x0);n=1;tol=1;whiletolepsx0=r;r=x0-myf(x0)/c;tol=norm(r-x0);n=n+1;if(n100000)disp('‘迭代步数太多，可能不收敛！');return;endend举例说明：解：首先建立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个多项式,542)(2341xxxxP,2)(2xxP,32)(23xxxP试进行下列操作：（1）求)(xP＝)(1xP＋)(2xP)(3xP（2）求)(xP的根（3）当x取矩阵A的每一个元素，求)(xP的值。其中：A＝5.2505.3275.04.12.11（4）当以矩阵A为自变量时，求)(xP的值。其中A的值与第（3）题相同。流程图源程序代码：%求积函数conv(x,y)%求平方根函数roots(p)%逐一取用矩阵中的数值函数polyval(p,x)%取用矩阵的函数polyvalm(p,A)是按照矩阵运算规则计算多项式的值%多项式p2，p3%矩阵A使用函数polyvalm(p,A)polyval(p,x)roots(p)conv(p2,p3)p1=[12405];p2=[12];p3=[123];x=[-11.2-1.4;0.7523.5;052.5];A=[-11.2-1.4;0.7523.5;052.5];开始结束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.1179x1=[0:0.05:1];y1=a(4)+a(3)*x1+a(2)*x1.^2+a(1)*x1.^3;plot(x,y,'*')holdonplot(x1,y1,'-r')运行结果：流程图：开始输入数据调用polyfit函数拟合令X1=[0:0.05:1]y1=a(4)+a(3)*x1+a(2)*x1.^2+a(1)*x1.^3;利用PLOT绘制图形结束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文件程序代码】functionjisuanx=[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)holdonsubplot(2,1,2)loglog(x0,y0)holdonjisuana=4.4717b=-0.9238流程图：开始输入数据调用polyfit函数拟合令A=exp(p(2))b=p(1)输入x0=0.0129:0.0001:2.46;y0=4.4717.*exp(-0.9238.*x0);利用PLOT绘制图形结束