matlab实验一:非线性方程求解-牛顿法
- 格式:doc
- 大小:22.00 KB
- 文档页数:2
非线性方程求根——牛顿迭代法一、牛顿迭代法的基本思想基本思想:将非线性方程逐步归结为某种线性方程求解。
设方程f (x )=0有近似根x k (f `(x k )≠0),将f (x )在x k 展开:(ξ在x 和x k 之间)2()()()()()()2!k k k k f f x f x f x x x x x ξ'''=+-+-()()()()k k k f x f x f x x x '≈+-可设记该线性方程的根为x k +1,则()()()0k k k f x f x x x '+-=1()()k k k k f x x x f x +=-'故f (x )=0可近似表示为即为Newton 法迭代格式。
(k =0,1,……)例:用Newton 迭代法求方程310x x --=在x 0=1.5附近的近似实根。
解:32()1,()31f x x x f x x '=--=-迭代公式为312131kk k k k x x x x x +--=--计算步骤如下:(1)取初值x 0=1.5;(2)按照迭代公式计算x 1;(3)若|x 1-x 0|<=0.00001,终止迭代;否则,x 0=x 1;转(2);(4)输出迭代次数和近似根.二、牛顿迭代法的实现MATLAB求解程序设计:方程及一阶导数函数:function[fun,dfun]=fun0(x)fun=x^3-x-1;%求原函数的值dfun=3*x^2-1;%求一阶导数的值计算主程序:clearx0=1.5;[fun,dfun]=fun0(x0);x1=x0-fun/dfun;i=1;while abs(x1-x0)>1e-5x0=x1;[fun,dfun]=fun0(x0);x1=x0-fun/dfun;i=i+1;enddisp('the solution is x1=')x1disp('the iter time is ')i计算结果为:the solution is x1=x1 =1.3247the iter time isi =4可见经过4次迭代即到达要求的精度,原方程的一个近似实数根为1.3247.三、牛顿迭代法的收敛性牛顿迭代法的迭代函数:)()()(x f x f x x '-=ϕ222)]([)()()]([)()()]([1)(x f x f x f x f x f x f x f x '''='''-'-='ϕ设f (x *)=0,f `(x *)≠0,则ϕ`(x *)=0,故Newton 迭代法在x *附近至少平方收敛。
一、概述MATLAB是一种强大的数学软件工具,它提供了许多优秀的数值计算和数据分析功能。
其中,牛顿拉夫逊法和快速分解法是两种常用的数值计算方法,它们在解决非线性方程组和矩阵分解等问题中具有重要的应用价值。
本文将介绍如何在MATLAB中实现这两种方法,并对它们的优缺点进行详细分析。
二、牛顿拉夫逊法的实现1. 算法原理牛顿拉夫逊法是一种用于求解非线性方程组的迭代算法。
它利用函数的一阶和二阶导数信息来不断逼近方程组的解,直到满足精度要求为止。
算法原理可以用以下公式表示:公式1其中,x表示解向量,F(x)表示方程组的函数向量,J(x)表示方程组的雅可比矩阵,δx表示解的更新量。
通过不断迭代更新x,最终得到方程组的解。
2. MATLAB代码实现在MATLAB中,可以通过编写函数来实现牛顿拉夫逊法。
以下是一个简单的示例代码:在这段代码中,首先定义了方程组的函数向量和雅可比矩阵,然后利用牛顿拉夫逊法进行迭代更新,直到满足精度要求为止。
通过这种方式,就可以在MATLAB中实现牛顿拉夫逊法,并应用于各种实际问题。
三、快速分解法的实现1. 算法原理快速分解法是一种用于矩阵分解的高效算法。
它利用矩阵的特定性质,通过分解为更小的子问题来加速计算过程。
算法原理可以用以下公式表示:公式2其中,A表示要分解的矩阵,L和U分别表示矩阵的下三角和上三角分解。
通过这种分解方式,可以将原始矩阵的计算量大大减小,提高求解效率。
2. MATLAB代码实现在MATLAB中,可以利用内置函数来实现快速分解法。
以下是一个简单的示例代码:在这段代码中,利用MATLAB内置的lu函数进行LU分解,得到矩阵的下三角和上三角分解。
通过这种方式,就可以在MATLAB中实现快速分解法,并应用于各种矩阵计算问题。
四、方法比较与分析1. 算法复杂度牛顿拉夫逊法和快速分解法在计算复杂度上有所不同。
牛顿拉夫逊法的迭代次数取决于所求解问题的非线性程度,通常需要较多的迭代次数。
⽜顿迭代法解⾮线性⽅程组(MATLAB版)⽜顿迭代法,⼜名切线法,这⾥不详细介绍,简单说明每⼀次⽜顿迭代的运算:⾸先将各个⽅程式在⼀个根的估计值处线性化(泰勒展开式忽略⾼阶余项),然后求解线性化后的⽅程组,最后再更新根的估计值。
下⾯以求解最简单的⾮线性⼆元⽅程组为例(平⾯⼆维定位最基本原理),贴出源代码:1、新建函数fun.m,定义⽅程组1 function f=fun(x);2 %定义⾮线性⽅程组如下3 %变量x1 x24 %函数f1 f25 syms x1 x26 f1 = sqrt((x1-4)^2 + x2^2)-sqrt(17);7 f2 = sqrt(x1^2 + (x2-4)^2)-5;8 f=[f1 f2];2、新建dfun.m,求出⼀阶微分⽅程1 function df=dfun(x);2 f=fun(x);3 df=[diff(f,'x1');diff(f,'x2')]; %雅克⽐矩阵3、建⽴newton.m,执⾏⽜顿迭代过程1 clear;clc2 format;3 x0=[0 0]; % 迭代初始值4 eps = 0.00001; % 定位精度要求5for i = 1:106 f = double(subs(fun(x0),{'x1''x2'},{x0(1) x0(2)}));7 df = double(subs(dfun(x0),{'x1''x2'},{x0(1) x0(2)})); % 得到雅克⽐矩阵8 x = x0 - f/df;9if(abs(x-x0) < eps)10break;11 end12 x0 = x; % 更新迭代结果13 end14 disp('定位坐标:');15 x16 disp('迭代次数:');17 i结果如下:定位坐标:x =0.0000 -1.0000迭代次数:i =4。
第7章 求解非线性方程7.1 多项式运算在MATLAB 中的实现一、多项式的表达n 次多项式表达为:n a +⋯⋯++=x a x a x a p(x )1-n 1-n 1n 0,是n+1项之和在MATLAB 中,n 次多项式可以用n 次多项式系数构成的长度为n+1的行向量表示[a0, a1,……an-1,an]二、多项式的加减运算设有两个多项式n a +⋯⋯++=x a x a x a p1(x )1-n 1-n 1n 0和m b +⋯⋯++=x b x b x b p2(x )1-m 1-m 1m 0。
它们的加减运算实际上就是它们的对应系数的加减运算。
当它们的次数相同时,可以直接对多项式的系数向量进行加减运算。
当它们的次数不同时,应该把次数低的多项式无高次项部分用0系数表示。
例2 计算()()1635223-+++-x x x xa=[1, -2, 5, 3]; b=[0, 0, 6, -1]; c=a+b例 3 设()6572532345++-+-=x x x x x x f ,()3532-+=x x x g ,求f(x)+g(x)f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3]; g1=[0, 0, 0, g];%为了和f 的次数找齐 f+g1, f-g1三、多项式的乘法运算conv(p1,p2)例4 在上例中,求f(x)*g(x)f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3];conv(f, g)四、多项式的除法运算[Q, r]=deconv(p1, p2)表示p1除以p2,给出商式Q(x),余式r(x)。
Q,和r 仍为多项式系数向量 例4 在上例中,求f(x)/g(x)f=[3, -5, 2, -7, 5, 6]; g=[3, 5, -3];[Q, r]=deconv(f, g)五、多项式的导函数p=polyder(P):求多项式P 的导函数p=polyder(P ,Q):求P ·Q 的导函数[p,q]=polyder(P ,Q):求P/Q 的导函数,导函数的分子存入p ,分母存入q 。
一、引言在数值计算中,求解非线性方程是一项常见的任务。
牛顿迭代法是一种常用且有效的方法,它通过不断逼近函数的零点来求解方程。
而在MATLAB中,我们可以利用其强大的数值计算功能来实现牛顿迭代法,快速求解各种非线性方程。
二、牛顿迭代法原理与公式推导1. 牛顿迭代法原理牛顿迭代法是一种利用函数的导数信息不断逼近零点的方法。
其核心思想是利用当前点的切线与x轴的交点来更新下一次迭代的值,直至逼近方程的根。
2. 公式推导与迭代过程假设要求解方程f(x)=0,在初始值x0附近进行迭代。
根据泰勒展开,对f(x)进行一阶泰勒展开可得:f(x) ≈ f(x0) + f'(x0)(x - x0)令f(x)≈0,则有:x = x0 - f(x0)/f'(x0)将x带入f(x)的表达式中,即得到下一次迭代的值x1:x1 = x0 - f(x0)/f'(x0)重复以上过程,直至达到精度要求或者迭代次数上限。
三、MATLAB中的牛顿迭代法实现1. 编写函数在MATLAB中,我们可以编写一个函数来实现牛顿迭代法。
需要定义原方程f(x)的表达式,然后计算其一阶导数f'(x)的表达式。
按照上述推导的迭代公式,编写循环语句进行迭代计算,直至满足精度要求或者达到最大迭代次数。
2. 调用函数求解方程在编写好牛顿迭代法的函数之后,可以通过在MATLAB命令窗口中调用该函数来求解具体的方程。
传入初始值、精度要求和最大迭代次数等参数,即可得到方程的近似根。
四、牛顿迭代法在工程实践中的应用1. 求解非线性方程在工程领域,很多问题都可以转化为非线性方程的求解问题,比如电路分析、控制系统设计等。
利用牛顿迭代法可以高效地求解这些复杂方程,为工程实践提供了重要的数值计算手段。
2. 优化问题的求解除了求解非线性方程外,牛顿迭代法还可以应用于优化问题的求解。
通过求解目标函数的导数等于0的方程,可以找到函数的极值点,从而解决各种优化问题。
非线性方程组求解的牛顿迭代法用MATLAB实现首先,我们需要定义非线性方程组。
假设我们要求解方程组:```f1(x1,x2)=0f2(x1,x2)=0```其中,`x1`和`x2`是未知数,`f1`和`f2`是非线性函数。
我们可以将这个方程组表示为向量的形式:```F(x)=[f1(x1,x2);f2(x1,x2)]=[0;0]```其中,`F(x)`是一个列向量。
为了实现牛顿迭代法,我们需要计算方程组的雅可比矩阵。
雅可比矩阵是由方程组的偏导数组成的矩阵。
对于方程组中的每个函数,我们可以计算其对每个变量的偏导数,然后将这些偏导数组成一个矩阵。
在MATLAB中,我们可以使用`jacobi`函数来计算雅可比矩阵。
以下是一个示例函数的定义:```matlabfunction J = jacobi(x)x1=x(1);x2=x(2);J = [df1_dx1, df1_dx2; df2_dx1, df2_dx2];end```其中,`x`是一个包含未知数的向量,`df1_dx1`和`df1_dx2`是`f1`对`x1`和`x2`的偏导数,`df2_dx1`和`df2_dx2`是`f2`对`x1`和`x2`的偏导数。
下一步是实现牛顿迭代法。
牛顿迭代法的迭代公式为:```x(k+1)=x(k)-J(x(k))\F(x(k))```其中,`x(k)`是第`k`次迭代的近似解,`\`表示矩阵的求逆操作。
在MATLAB中,我们可以使用如下代码来实现牛顿迭代法:```matlabfunction x = newton_method(x_initial)max_iter = 100; % 最大迭代次数tol = 1e-6; % 收敛阈值x = x_initial; % 初始解for k = 1:max_iterF=[f1(x(1),x(2));f2(x(1),x(2))];%计算F(x)J = jacobi(x); % 计算雅可比矩阵 J(x)delta_x = J \ -F; % 计算增量 delta_xx = x + delta_x; % 更新 xif norm(delta_x) < tolbreak; % 达到收敛条件,停止迭代endendend```其中,`x_initial`是初始解的向量,`max_iter`是最大迭代次数,`tol`是收敛阈值。
matlab 实现牛顿迭代法求解非线性方程组已知非线性方程组如下3*x1-cos(x2*x3)-1/2=0x1A2-81*(x2+A2+si n(x3)+=0exp(-x1*x2)+20*x3+(10*pi-3)/3=0求解要求精度达到首先建立函数fun 储存方程组编程如下将保存到工作路径中function f=fun(x);%定义非线性方程组如下%变量x1 x2 x3%函数f1 f2 f3syms x1 x2 x3 f1=3*x1-cos(x2*x3)-1/2;f2=x1A2-81*(x2+A2+sin(x3)+; f3=exp(-x1*x2)+20*x3+(10*pi-3)/3;f=[f1 f2 f3];建立函数dfun 用来求方程组的雅克比矩阵将保存到工作路径中function df=dfun(x); %用来求解方程组的雅克比矩阵储存在dfun 中f=fun(x);df=[diff(f,'x1');diff(f,'x2');diff(f,'x3')];df=conj(df');编程牛顿法求解非线性方程组将保存到工作路径中function x=newton(x0,eps,N);con=0;eps 为精度要求 N 为最大迭代步数 con 用来记录结果是否收敛 for i=1:N;x=x0-f/df;for j=1: length(x0);il(i,j)=x(j);end if norm(x-x0)<eps con=1;break;end x0=x;endfid=fopen('','w');fprintf(fid,'iteration');for j=1:length(x0)endfor j=1:ifor k=1:length(x0) %其中 x0 为迭代初值 f=subs(fun(x0),{'x1' 'x2' 'x3'},{x0(1)x0(2) x0(3)}); df=subs(dfun(x0),{'x1' 'x2' 'x3'},{x0(1)x0(2) x0(3)}); %以下是将迭代过程写入 txt 文档文件名为fprintf(fid,' x%d',j);fprintf(fid,'\n%6d ',j);fprintf(fid,' %',il(j,k));end endif con==1fprin tf(fid,'\n 计算结果收敛!');endif con==0fprin tf(fid,'\n 迭代步数过多可能不收敛!');endfclose(fid);运行程序在matlab中输入以下内容newt on([ ],,20)输出结果»nawtnnC[0. 1 0. 1 -0, 11,0.00001,20)ans =0. 5000 0. OOCO -0. 523S在iteration 中查看迭代过程iteration x1 x2 x31 2j.1:erat ion xl X30. 031238-0.51960110.490718320. 5090110. 003498-0.521634 430. 5009 2&0. 000756-0.523391 540. 5002E70. 000076-0.52355065)5000190+ 000013-0. 5235947S0. 50LI005 D. 000002-0.S235983r CL 5000000+OOOCOO-0. 5235999 计尊结果收敘!.mulStablePoint 用不动点迭代法求非线性方程组的一个根function [r,n]=mulStablePoint(F,x0,eps)%非线性方程组:f%初始解:a%解的精度:eps%求得的一组解:r%迭代步数:nif nargin==2eps=;endx0 = transpose(x0);n=1;tol=1;while tol>epsr= subs(F,findsym(F),x0);%迭代公式tol=norm(r-x0); %注意矩阵的误差求法,norm 为矩阵的欧几里德范数n=n+1;x0=r;if(n>100000)%迭代步数控disp('迭代步数太多,可能不收敛!');return;endendxO=[O 0 0];[r,n, data]=budo ng(xO);disp('不动点计算结果为')x仁[1 1 1];x2=[2 2 2];[x, n,data]=n ew_ton (x0);disp('初始值为0,牛顿法计算结果为:’)[x, n,data]=n ew_ton (x1);disp('初始值为1,牛顿法计算结果为:')[x, n,data]=n ew_ton (x2);disp ('初始值为2,牛顿法计算结果为:')function[r,n, data]=bud on g(x0, tol) if nargin=-1tol=1e-3 :endx1=budo ng fun( x0) ;n=1;while( norm(x1-x0))tol)&(n500)xO=x1;x1=budo ng_fun( x0);n=n+1:data(:, n)=x1 ;endr=x1 :function [x,n,data]=new_ton(x0, tol)if nargin=-1tol=1e-8 ;endx1=x0-budo ng_fun (x0)/df1(x0);n=1;while (n orm(x1-x0))tol)x0=x1;x1=x0-budo ng_fun (x0)/df1(x0);n=n+1;data(:, n)=x1;endx=x1;function f=bud ong_fun(x)f(1)=3* x(1)-cos(x(2)*x (3))-1/2;f(2)=x(1)A2-81*(x(2)+A2+si n(x(3))+; f(3)=exp(-x(1)*x(2))+20* x(3)+10* pi/3-1; f=[f(1)*f(2)*f(3)] ;fun ctio n f=df1(x)f=[3s in (x(2)*x(3))*x(3) si n(x(2)*x(3))*x(2)2* x(1)-162*(x(2)+cos(x(3))exp(-x(1)*x(2))*(-x (2))exp(-x(1)*x(2))*(-x(1))20]; 结果:不动点计算结果为r=+012*NaN -Inf初始值为0,牛顿法计算结果为:x=初始值为1,牛顿法计算结果为:x=初始值为2,牛顿法计算结杲为:x=。
一、简介Matlab是一种十分常用的科学计算软件,其功能强大,可以进行各种数值计算、数据分析和可视化操作。
而牛顿迭代法是一种用于求解方程的数值算法,可以有效地计算出函数的根。
本文将重点介绍如何使用Matlab进行牛顿迭代法来计算重根。
二、牛顿迭代法原理1. 牛顿迭代法是一种迭代逼近的方法,通过不断迭代得到更接近函数零点的近似值。
其公式如下:X_{n+1} = X_n - \frac{f(X_n)}{f'(X_n)}其中,X_{n+1}为下一次迭代的近似值,X_n为当前的近似值,f(X)为函数值,f'(X)为函数的导数值。
2. 牛顿迭代法的优点是收敛速度快,而缺点是对初始值的选择敏感,可能会产生不收敛的情况。
三、在Matlab中使用牛顿迭代法1. 在Matlab中,可以使用内置的函数`fzero`来进行牛顿迭代法的计算。
其语法如下:x = fzero(fun,x0)其中,fun为要求解的函数句柄,x0为起始点的初始值,x为函数的根。
2. 需要注意的是,在使用`fzero`函数时,需要提供函数的句柄,即在Matlab中定义要求解的函数,并使用`(x)`符号来表示函数的自变量。
另外,还需要提供初始值x0,可以根据具体问题来选择较为合适的初始值。
3. 以下是一个简单的使用牛顿迭代法求解函数根的示例代码:```matlabf = (x) x^3 - 2*x - 5;x0 = 2;x = fzero(f, x0);disp(x);```四、示例接下来,我们将通过一个具体的示例来演示如何使用Matlab的牛顿迭代法来计算重根。
1. 问题描述假设有如下方程:f(x) = x^3 - 2x^2 + 3x - 6我们希望使用牛顿迭代法来计算函数f(x)的重根。
2. 解决过程在Matlab中定义函数f(x):```matlabf = (x) x^3 - 2*x^2 + 3*x - 6;```选择初始值x0,并利用`fzero`函数进行牛顿迭代法的计算:```matlabx0 = 2;x = fzero(f, x0);disp(x);```3. 结果分析经过计算,可以得到函数f(x)的一个重根为x=2.这样,我们就成功地使用Matlab的牛顿迭代法来计算重根。
非线性方程组求解1.mulStablePoint用不动点迭代法求非线性方程组的一个根function [r,n]=mulStablePoint(F,x0,eps)%非线性方程组:f%初始解:a%解的精度:eps%求得的一组解:r%迭代步数:nif nargin==2eps=1.0e-6;endx0 = transpose(x0);n=1;tol=1;while tol>epsr= subs(F,findsym(F),x0); %迭代公式tol=norm(r-x0); %注意矩阵的误差求法,norm为矩阵的欧几里德范数n=n+1;x0=r;if(n>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endend2.mulNewton用牛顿法法求非线性方程组的一个根function [r,n]=mulNewton(F,x0,eps)if nargin==2eps=1.0e-4;endx0 = transpose(x0);Fx = subs(F,findsym(F),x0);var = findsym(F);dF = Jacobian(F,var);dFx = subs(dF,findsym(dF),x0);r=x0-inv(dFx)*Fx;n=1;tol=1;while tol>epsx0=r;Fx = subs(F,findsym(F),x0);dFx = subs(dF,findsym(dF),x0);r=x0-inv(dFx)*Fx; %核心迭代公式tol=norm(r-x0);n=n+1;if(n>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endend3.mulDiscNewton用离散牛顿法法求非线性方程组的一个根function [r,m]=mulDiscNewton(F,x0,h,eps)format long;if nargin==3eps=1.0e-8;endn = length(x0);fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endr=transpose(x0)-inv(J)*fx;m=1;tol=1;while tol>epsxs=r;fx = subs(F,findsym(F),xs);J = zeros(n,n);for i=1:nx1 = xs;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endr=xs-inv(J)*fx; %核心迭代公式tol=norm(r-xs);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;4.mulMix用牛顿-雅可比迭代法求非线性方程组的一个根function [r,m]=mulMix(F,x0,h,l,eps)if nargin==4eps=1.0e-4;endn = length(x0);J = zeros(n,n);Fx = subs(F,findsym(F),x0);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));C =D - J;inD = inv(D);H = inD*C;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = Hm*inD*Fx;r = transpose(x0)-dr; m=1;tol=1;while tol>epsx0=r;Fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));C =D - J;inD = inv(D);H = inD*C;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = Hm*inD*Fx;r = x0-dr; %核心迭代公式tol=norm(r-x0);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endend5.mulNewtonSOR用牛顿-SOR迭代法求非线性方程组的一个根function [r,m]=mulNewtonSOR(F,x0,w,h,l,eps)if nargin==5eps=1.0e-4;endn = length(x0);J = zeros(n,n);Fx = subs(F,findsym(F),x0);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));L = -tril(J-D);U = -triu(J-D);inD = inv(D-w*L);H = inD*(D - w*D+w*L);;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = w*Hm*inD*Fx;r = transpose(x0)-dr;m=1;tol=1;while tol>epsx0=r;Fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endD = diag(diag(J));L = -tril(J-D);U = -triu(J-D);inD = inv(D-w*L);H = inD*(D - w*D+w*L);;Hm = eye(n,n);for i=1:l-1Hm = Hm + power(H,i);enddr = w*Hm*inD*Fx;r = x0-dr; %核心迭代公式tol=norm(r-x0);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endend6.mulDNewton用牛顿下山法求非线性方程组的一个根function [r,m]=mulDNewton(F,x0,eps)%非线性方程组:F%初始解:x0%解的精度:eps%求得的一组解:r%迭代步数:nif nargin==2eps=1.0e-4;endx0 = transpose(x0);dF = Jacobian(F);m=1;tol=1;while tol>epsttol=1;w=1;Fx = subs(F,findsym(F),x0);dFx = subs(dF,findsym(dF),x0);F1=norm(Fx);while ttol>=0 %下面的循环是选取下山因子w的过程r=x0-w*inv(dFx)*Fx; %核心的迭代公式Fr = subs(F,findsym(F),r);ttol=norm(Fr)-F1;w=w/2;endtol=norm(r-x0);m=m+1;x0=r;if(m>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endend7.mulGXF1用两点割线法的第一种形式求非线性方程组的一个根function [r,m]=mulGXF1(F,x0,x1,eps)format long;if nargin==3eps=1.0e-4;endx0 = transpose(x0);x1 = transpose(x1);n = length(x0);fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);for i=1:nxt = x1;xt(i) = x0(i);J(:,i) = (subs(F,findsym(F),xt)-fx1)/h(i);endr=x1-inv(J)*fx1;m=1;tol=1;while tol>epsx0 = x1;x1 = r;fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);for i=1:nxt = x1;xt(i) = x0(i);J(:,i) = (subs(F,findsym(F),xt)-fx1)/h(i);endr=x1-inv(J)*fx1;tol=norm(r-x1);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;8.mulGXF2用两点割线法的第二种形式求非线性方程组的一个根function [r,m]=mulGXF2(F,x0,x1,eps)format long;if nargin==3eps=1.0e-4;endx0 = transpose(x0);x1 = transpose(x1);n = length(x0);fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);xt = x1;xt(1) = x0(1);J(:,1) = (subs(F,findsym(F),xt)-subs(F,findsym(F),x1))/h(1);for i=2:nxt = x1;xt(1:i) = x0(1:i);xt_m = x1;xt_m(1:i-1) = x0(1:i-1);J(:,i) = (subs(F,findsym(F),xt)-subs(F,findsym(F),xt_m))/h(i);endr=x1-inv(J)*fx1;m=1;tol=1;while tol>epsx0 = x1;x1 = r;fx = subs(F,findsym(F),x0);fx1 = subs(F,findsym(F),x1);h = x0 - x1;J = zeros(n,n);xt = x1;xt(1) = x0(1);J(:,1) = (subs(F,findsym(F),xt)-subs(F,findsym(F),x1))/h(1);for i=2:nxt = x1;xt(1:i) = x0(1:i);xt_m = x1;xt_m(1:i-1) = x0(1:i-1);J(:,i) = (subs(F,findsym(F),xt)-subs(F,findsym(F),xt_m))/h(i);endr=x1-inv(J)*fx1;tol=norm(r-x1);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;9.mulVNewton用拟牛顿法求非线性方程组的一组解function [r,m]=mulVNewton(F,x0,A,eps)%方程组:F%方程组的初始解:x0% 初始A矩阵:A%解的精度:eps%求得的一组解:r%迭代步数:mif nargin==2A=eye(length(x0)); %A取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendx0 = transpose(x0);Fx = subs(F, findsym(F),x0);r=x0-A\Fx;m=1;tol=1;while tol>epsx0=r;Fx = subs(F, findsym(F),x0);r=x0-A\Fx;y=r-x0;Fr = subs(F, findsym(F),r);z= Fr-Fx;A1=A+(z-A*y)*transpose(y)/norm(y); %调整A A=A1;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endtol=norm(r-x0);end10.mulRank1用对称秩1算法求非线性方程组的一个根function [r,n]=mulRank1(F,x0,A,eps)if nargin==2l = length(x0);A=eye(l); %A取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendfx = subs(F,findsym(F),x0);r=transpose(x0)-inv(A)*fx;n=1;tol=1;while tol>epsx0=r;fx = subs(F,findsym(F),x0);r=x0-inv(A)*fx;y=r-x0;fr = subs(F,findsym(F),r);z = fr-fx;A1=A+ fr *transpose(fr)/(transpose(fr)*y); %调整A A=A1;n=n+1;if(n>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endtol=norm(r-x0);end11.mulDFP用D-F-P算法求非线性方程组的一组解function [r,n]=mulDFP(F,x0,A,eps)if nargin==2l = length(x0);B=eye(l); %A取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendfx = subs(F,findsym(F),x0);r=transpose(x0)-B*fx;n=1;tol=1;while tol>epsx0=r;fx = subs(F,findsym(F),x0);r=x0-B*fx;y=r-x0;fr = subs(F,findsym(F),r);z = fr-fx;B1=B+ y*y'/(y'*z)-B*z*z'*B/(z'*B*z); %调整AB=B1;n=n+1;if(n>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endtol=norm(r-x0);end12.mulBFS用B-F-S算法求非线性方程组的一个根function [r,n]=mulBFS(F,x0,B,eps)if nargin==2l = length(x0);B=eye(l); %B取为单位阵eps=1.0e-4;elseif nargin==3eps=1.0e-4;endendfx = subs(F,findsym(F),x0);r=transpose(x0)-B*fx;n=1;tol=1;while tol>epsx0=r;fx = subs(F,findsym(F),x0);r=x0-B*fx;y=r-x0;fr = subs(F,findsym(F),r);z = fr-fx;u = 1 + z'*B*z/(y'*z);B1= B+ (u*y*y'-B*z*y'-y*z'*B)/(y'*z); %调整B B=B1;n=n+1;if(n>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endtol=norm(r-x0);end13.mulNumYT用数值延拓法求非线性方程组的一组解function [r,m]=mulNumYT(F,x0,h,N,eps)format long;if nargin==4eps=1.0e-8;endn = length(x0);fx0 = subs(F,findsym(F),x0);x0 = transpose(x0);J = zeros(n,n);for k=0:N-1fx = subs(F,findsym(F),x0);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endinJ = inv(J);r=x0-inJ*(fx-(1-k/N)*fx0);x0 = r;endm=1;tol=1;while tol>epsxs=r;fx = subs(F,findsym(F),xs);J = zeros(n,n);for i=1:nx1 = xs;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-fx)/h(i);endr=xs-inv(J)*fx; %核心迭代公式tol=norm(r-xs);m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;14.DiffParam1用参数微分法中的欧拉法求非线性方程组的一组解function r=DiffParam1(F,x0,h,N)%非线性方程组:f%初始解:x0%数值微分增量步大小:h%雅可比迭代参量:l%解的精度:eps%求得的一组解:r%迭代步数:nx0 = transpose(x0);n = length(x0);ht = 1/N;Fx0 = subs(F,findsym(F),x0);for k=1:NFx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h(i);J(:,i) = (subs(F,findsym(F),x1)-Fx)/h(i);endinJ = inv(J);r = x0 - ht*inJ*Fx0;x0 = r;end15.DiffParam2用参数微分法中的中点积分法求非线性方程组的一组解function r=DiffParam2(F,x0,h,N)%非线性方程组:f%初始解:x0%数值微分增量步大小:h%雅可比迭代参量:l%解的精度:eps%求得的一组解:r%迭代步数:nx0 = transpose(x0);n = length(x0);ht = 1/N;Fx0 = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nxt = x0;xt(i) = xt(i)+h(i);J(:,i) = (subs(F,findsym(F),xt)-Fx0)/h(i);endinJ = inv(J);x1 = x0 - ht*inJ*Fx0;for k=1:Nx2 = x1 + (x1-x0)/2;Fx2 = subs(F,findsym(F),x2);J = zeros(n,n);for i=1:nxt = x2;xt(i) = xt(i)+h(i);J(:,i) = (subs(F,findsym(F),xt)-Fx2)/h(i);endinJ = inv(J);r = x1 - ht*inJ*Fx0;x0 = x1;x1 = r;end16.mulFastDown用最速下降法求非线性方程组的一组解function [r,m]=mulFastDown(F,x0,h,eps)format long;if nargin==3eps=1.0e-8;endn = length(x0);x0 = transpose(x0);m=1;tol=1;while tol>epsfx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;J(:,i) = (subs(F,findsym(F),x1)-fx)/h;endlamda = fx/sum(diag(transpose(J)*J));r=x0-J*lamda; %核心迭代公式fr = subs(F,findsym(F),r);tol=dot(fr,fr);x0 = r;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;17.mulGSND用高斯牛顿法求非线性方程组的一组解function [r,m]=mulGSND(F,x0,h,eps)format long;if nargin==3eps=1.0e-8;endn = length(x0);x0 = transpose(x0);m=1;tol=1;while tol>epsfx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;J(:,i) = (subs(F,findsym(F),x1)-fx)/h;endDF = inv(transpose(J)*J)*transpose(J);r=x0-DF*fx; %核心迭代公式tol=norm(r-x0);x0 = r;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;18.mulConj用共轭梯度法求非线性方程组的一组解function [r,m]=mulConj(F,x0,h,eps)format long;if nargin==3eps=1.0e-6;endn = length(x0);x0 = transpose(x0);fx0 = subs(F,findsym(F),x0);p0 = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)*(1+h);p0(:,i) = -(subs(F,findsym(F),x1)-fx0)/h;endm=1;tol=1;while tol>epsfx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;J(:,i) = (subs(F,findsym(F),x1)-fx)/h;endlamda = fx/sum(diag(transpose(J)*J));r=x0+p0*lamda; %核心迭代公式fr = subs(F,findsym(F),r);Jnext = zeros(n,n);for i=1:nx1 = r;x1(i) = x1(i)+h;Jnext(:,i) = (subs(F,findsym(F),x1)-fr)/h;endabs1 = transpose(Jnext)*Jnext;abs2 = transpose(J)*J;v = abs1/abs2;if (abs(det(v)) < 1)p1 = -Jnext+p0*v;elsep1 = -Jnext;endtol=norm(r-x0);p0 = p1;x0 = r;m=m+1;if(m>100000) %迭代步数控制 disp('迭代步数太多,可能不收敛!');return;endendformat short;19.mulDamp用阻尼最小二乘法求非线性方程组的一组解function [r,m]=mulDamp(F,x0,h,u,v,eps)format long;if nargin==5eps=1.0e-6;endFI = transpose(F)*F/2;n = length(x0);x0 = transpose(x0);m=1;tol=1;while tol>epsj = 0;fx = subs(F,findsym(F),x0);J = zeros(n,n);for i=1:nx1 = x0;x1(i) = x1(i)+h;afx = subs(F,findsym(F),x1);J(:,i) = (afx-fx)/h;endFIx = subs(FI,findsym(FI),x0);for i=1:nx2 = x0;x2(i) = x2(i)+h;gradFI(i,1) = (subs(FI,findsym(FI),x2)-FIx)/h;ends=0;while s==0A = transpose(J)*J+u*eye(n,n);p = -A\gradFI;r = x0 + p;FIr = subs(FI,findsym(FI),r);if FIr<FIxif j == 0u = u/v;j = 1;elses=1;endelseu = u*v;j = 1;if norm(r-x0)<epss=1;endendendx0 = r;tol = norm(p);m=m+1;if(m>100000) %迭代步数控制disp('迭代步数太多,可能不收敛!');return;endendformat short;。
实验一:非线性方程求解
程序1:二分法:
syms f x;
f=input('请输入f(x)=');
A=input('请输入根的估计范围[a,b]='); e=input('请输入根的误差限e='); while (A(2)-A(1))>e
c=(A(1)+A(2))/2;
x=A(1);
f1=eval(f);
x=c;
f2=eval(f);
if (f1*f2)>0
A(1)=c;
else
A(2)=c;
end
end
c=(A(1)+A(2))/2;
fprintf('c=%.6f\na=%.6f\nb=%.6f\n',c,A)
用二分法计算方程:
1.请输入f(x)=sin(x)-x^2/2
请输入根的估计范围[a,b]=[1,2]
请输入根的误差限e=0.5e-005
c=1.404413
a=1.404411
b=1.404415
2.请输入f(x)=x^3-x-1
请输入根的估计范围[a,b]=[1,1.5]
请输入根的误差限e=0.5e-005
c=1.324717
a=1.324715
b=1.324718
程序2:newton法:
syms f x;
f=input('请输入f(x)=');
df=diff(f); x0=input('请输入迭代初值x0=');
e1=input('请输入奇异判断e1=');
e2=input('请输入根的误差限e2=');
N=input('请输入迭代次数限N=');
k=1;
while (k<N)
x=x0;
if abs(eval(f))<e1
fprintf('奇异!\nx=%.6f\n迭代次数为:%d\n',x0,k)
break
else
x1=x0-eval(f)/eval(df);
if abs(x1-x0)<e2
fprintf('x=%.6f\n迭代次数为:%d\n',x1,k)
break
else
x0=x1;
k=k+1;
end
end
end
if k>=N
fprintf('失败\n')
end
用newton法计算方程:
1.请输入f(x)=x*exp(x)-1
请输入迭代初值x0=0.5
请输入奇异判断e1=0.1e-010
请输入根的误差限e2=0.5e-005
请输入迭代次数限N=10
x=0.567143
迭代次数为:4
2.请输入f(x)=x^3-x-1
请输入迭代初值x0=1
请输入奇异判断e1=0.1e-010
请输入根的误差限e2=0.5e-005
请输入迭代次数限N=10
x=1.324718
迭代次数为:5
3.1:请输入f(x)=(x-1)^2*(2*x-1)
请输入迭代初值x0=0.45
请输入奇异判断e1=0.1e-010
请输入根的误差限e2=0.5e-005
请输入迭代次数限N=10
x=0.500000
迭代次数为:4
3.2:请输入f(x)=(x-1)^2*(2*x-1)
请输入迭代初值x0=0.65
请输入奇异判断e1=0.1e-010
请输入根的误差限e2=0.5e-005
请输入迭代次数限N=10
x=0.500000
迭代次数为:9
3.3:请输入f(x)=(x-1)^2*(2*x-1)
请输入迭代初值x0=0.55
请输入奇异判断e1=0.1e-010
请输入根的误差限e2=0.5e-005
请输入迭代次数限N=10
x=0.500000
迭代次数为:4
程序3:改进的newton法:
syms f x;
f=input('请输入f(x)=');
df=diff(f);
x0=input('请输入迭代初值x0=');
e1=input('请输入奇异判断e1=');
e2=input('请输入根的误差限e2=');
N=input('请输入迭代次数限N=');
k=1;
while (k<N)
x=x0;
if abs(eval(f))<e1
fprintf('奇异!\nx=%.6f\n迭代次数为:%d\n',x0,k)
break
else
x1=x0-2*eval(f)/eval(df);
if abs(x1-x0)<e2
fprintf('x=%.6f\n迭代次数为:%d\n',x1,k)
break
else
x0=x1;
k=k+1;
end
end
end
if k>=N
fprintf('失败\n')
end
用改进的newton法计算方程:
1.请输入f(x)=(x-1)^2*(2*x-1)
请输入迭代初值x0=0.55
请输入奇异判断e1=0.1e-010
请输入根的误差限e2=0.5e-005
请输入迭代次数限N=10
失败
2.请输入f(x)=(x-1)^2*(2*x-1)
请输入迭代初值x0=0.55
请输入奇异判断e1=0.1e-010
请输入根的误差限e2=0.5e-005
请输入迭代次数限N=20
失败
3.请输入f(x)=(x-1)^2*(2*x-1)
请输入迭代初值x0=0.55
请输入奇异判断e1=0.1e-010
请输入根的误差限e2=0.5e-005
请输入迭代次数限N=100
失败。