湖北民族学院《数值分析》实验报告
实验名称 解非线性方程的数值计算方法 实验时间 2011年 11月 9日 姓名
王亚敏
班级
0209408
学号
020940807
成绩
实验报告内容要求：
一、实验目的与要求；二、实验内容；三、算法描述(数学原理或设计思路、计算公式、计算步骤)； 四、程序代码；五、数值结果；六、计算结果分析（如初值对结果的影响；不同方法的比较；该方法的特点和改进等）；七、实验中出现的问题，解决方法及体会（整个实验过程中（包括程序编写，上机调试等）出现的问题及其处理等广泛的问题）.
一、实验目的：
① 掌握二分法、不动点迭代法、牛顿迭代法、Steffensen 迭代法等常用的非线性方程迭代算法； ② 培养编程与上机调试能力. 二、实验内容：
1．请分别用二分法、不动点迭代法、牛顿迭代法、Steffensen 迭代法求方程
32
()330f x x x x =+--=在1.5 附近的根.
参考答案：原方程在1.5 附近的根为1.732051x =
2．请自选两个非线性方程，分别用二分法、不动点迭代法、牛顿迭代法、Steffensen 迭代法求解. 三、算法
（一） 1．二分法
计算()0f x =的二分法如下：
① 输入求根取间[,]a b 和误差控制量ε，定义函数. 如果 ()()0f a f b <，转②；否则退出选用其它求根方法 ② 当a b ε->时，计算中点()/2x a b =+以及的值； 分情况处理
()f x ε<：停止计算，，转④
()()0f a f x <：修正区间[,][,]ax ab → ()()0f x f b <：修正区间[,][,]xb ab →
③ *2a b
x +=
④ 输出近似根*
x
2．不动点(简单)迭代法算法思想
先将f(x)=0化成x=φ(x),然后取x 0,计算
x n+1=φ(x n ) n=0,1,2,…
若要求x *
满足f(x *
)=0,则x *
=)(*x ?，称x *
为函数)(x ?的一个不动点。
求 f(x)的零点就等价于求?（x ）的不动点，选择一个初始近似值x 0，将它代入x=?（x ）右 端，即可
求的)(01x x ?=可以如此反复迭代计算x n+1=φ(x n ),
?（x ）称为迭代函数，如此得到的序列{x n }有极限
*∞
→=x x k k lim 则称迭代方程收敛，且x *=)(*x ?为)(x ?的不动点，故称x n+1=φ(x n )为不动点迭代法。
3．牛顿迭代法
给定初始值x ，ε为根的容许误差，η为的容许误差，N 为迭代次数的容许值. ① 如果'0()0f x =或迭代次数大于N ，则算法失败，结束；否则执行② ② 计算010'0()
()
f x x x f x =-
③ 若10x x ε-<或1()f x η<，则输出x ，程序结束；否则执行④ ④ 令0
1x x =，转向①
4、Steffensen 迭代法
如果把Aitken 加速技巧与不动点迭代结合,则可得到如下的迭代:
k
k k k k k k k k k k x y z x y x x y z x y +---
===+2)(),
(),(21?? (k=0,1,2,........)
这称为Steffensen 迭代法.它是二阶收敛的.它让不收敛的收敛,即使是收敛的用Steffensen 后可达到二阶收敛.
四、程序代码
1、二分法
function x=erfenfa(f,a,b,n,delta) f=inline(f); %定义f. yb=feval(f,b);
disp('i x') %以指定格式输出'i','x' for i=1:n
ya=feval(f,a); %f 的在a 点的值给yb. s=(a+b)/2; %开始二分区间 c=feval(f,s); %计算 if c==0 %判断二分区间 break
elseif (c*ya)<0 b=s; else a=s; end
if abs(b-a) disp(sprintf('%d %10e',i,s)) %以以上规定格式输出解果 x=s; 2、迭代法 function y=diedaifa(f,x0,w,max) f=inline(f); disp('k y') disp(sprintf('%d %15.14f',0,x0)) for k=1:max y=feval(f,x0); if abs(y-x0) break end x0=y; disp(sprintf('%d %15.14f',k,y)) end k 3、牛顿迭代 function y=newton_c(f,x0,w,max) f=inline(f); disp('k y') disp(sprintf('%d %15.14f',0,x0)) for k=1:max y=feval(f,x0); if abs(y-x0) break end x0=y; disp(sprintf('%d %15.14f',k,y)) end k 4、Steffensen迭代法 function s=steff(f,x0,d,max) f=inline(f); %定义函数 x(1)=x0; %给初值 disp('k x y z'); %指定格式输出k,x,y,z for k=1:max y(k)=feval(f,x(k)); %本函数在它指定的点的值 z(k)=feval(f,y(k)); %本函数在它指定的点的值 x(k+1)=x(k)-(y(k)-x(k))^2/(z(k)-2*y(k)+x(k)); %算法公式 if abs(x(k+1)-x(k)) break disp(sprintf('%d %f %f %f',k,x(k),y(k),z(k))); %输出结果end s=x(k+1); %输出结果 五、计算结果 （一） 1、二分法 x=erfenfa('x^3+x^2-3*x-3',1.0,2.0,6,10^(-3)) i x 1 1.500000e+000 2 1.750000e+000 3 1.625000e+000 4 1.687500e+000 5 1.718750e+000 6 1.734375e+000 x = 1.7344 2、迭代法 y=diedaifa('(-x^2+3*x+3)^(1/3)',1.5, 10^(-6),9) k y 0 1.50000000000000 1 1.73801332244322 2 1.73173933329688 3 1.73206685838641 4 1.73204997984982 k = 5 y = 1.7321 3、牛顿迭代 y=newton_c('x-(x^3+x^2-3*x-3)/(3*x^2+2*x-3)',1.5, 10^(-6),9) k y 0 1.50000000000000 1 1.77777777777778 2 1.73336066694000 3 1.73205192940947 4 1.73205080756970 k = 5 y = 1.7321 4、Steffensen迭代法 s=steff('(-x^2+3*x+3)^(1/3)',1.5,10^(-6),15) k x y z 1 1.500000 1.738013 1.731739 2 1.731900 1.732059 1.732050 s = 1.7321 （二）方程x^3-x^2-1 1、二分法 x=erfenfa('x^3-x^2-1',1.0,2.0,6,10^(-3)) i x 1 1.500000e+000 2 1.250000e+000 3 1.375000e+000 4 1.437500e+000 5 1.468750e+000 6 1.453125e+000 x = 1.4531 2、不动点迭代 y=diedaifa('(x^2+1)^(1/3)',1.5, 10^(-6),9) k y 0 1.50000000000000 1 1.48124803420369 2 1.47270572963939 3 1.46881731366450 4 1.46704797320060 5 1.46624301011475 6 1.46587682016881 7 1.46571024077587 8 1.46563446523851 9 1.46559999585332 k = 9 y = 1.4656 3、牛顿迭代 y=newton_c('x-(x^3-x^2-1)/(3*x^2-2*x)',1.5, 10^(-6),9) k y 0 1.50000000000000 1 1.46666666666667 2 1.46557239057239 3 1.46557123187807 k = 4 y = 1.4656 4、Steffensen迭代法 s=steff('(x^2+1)^(1/3)',1.5,10^(-6),15) k x y z 1 1.500000 1.481248 1.472706 2 1.465558 1.465565 1.465569 s = 1.4656 （三）方程想x^3+4x^2-10 1、二分法 x=erfenfa('x^3+4*x^2-10',1.0,2.0,6,10^(-3)) i x 1 1.500000e+000 2 1.250000e+000 3 1.375000e+000 4 1.312500e+000 5 1.343750e+000 6 1.359375e+000 x = 1.3594 2、不动点迭代 y=diedaifa('(1/2)*((10-x^3)^(1/2))',1.4, 10^(-6),9) k y 0 1.40000000000000 1 1.34684817258665 2 1.37448330427545 3 1.36045214519208 4 1.36766537414766 5 1.36398040574705 6 1.36586902918221 7 1.36490266970374 8 1.36539755025157 9 1.36514422782316 k = 9 y = 1.3651 3、牛顿迭代 y=newton_c('x-(x^3+4*(x^2)-10)/(3*x^2+8*x)',1.5, 10^(-6),9) k y 0 1.50000000000000 1 1.37333333333333 2 1.36526201487463 3 1.36523001391615 k = 4 y = 1.3652 4、Steffensen迭代法 s=steff('(1/2)*((10-x^3)^(1/2))',1.5,10^(-6),15) k x y z 1 1.500000 1.286954 1.402541 2 1.361886 1.366937 1.364355 3 1.365228 1.365231 1.365230 s = 1.3652 六、计算结果分析 从计算结果可以看出，不动点迭代法、牛顿迭代法、Steffensen迭代法的结果都比二分法要精确。Steffensen迭代法的收敛速度是最快的，最慢的是不动点迭代法，其收敛速度由快到慢分别是Steffensen迭代法、牛顿迭代法、二分法、不动点迭代法。从整体上看，Steffensen迭代法的方法最快有比较精确，Steffensen 迭代法相对其他方法是最好的方法。 七、实验中出现的问题，解决方法及体会 若取得迭代公式不收敛导致计算结果出现很大的偏差。例如32 +-=,其迭代方程如果为 4100 x x 32 x x =-|’>1不收敛，其算得的结果是 104 =-在 1.5附近区间[1.4,1.6]之间|32 x x 104 y=diedaifa('((10-4*x^2)^(1/3))',1.4, 10^(-6),9) k y 0 1.40000000000000 1 1.29266081401913 2 1.49122462562746 3 1.03384070955176 4 1.78889183191316 5 0.70477483421324 6 2.46841603863811 7 1.59288056868044 8 2.88354123868011 9 2.28818326392814 k =9 y = 2.2882 + 2.1247i 它的值明显明显偏差很大，所以在选择方程式一定要选好恰当的。其次，就算选好方程，但是区间范围也 要把握好，例如32 10x x --=选择的迭代公式2 1 1x x =+ 在[1.4,1.5]收敛，但是在区间[0,3]之间就不收敛。 教 师 评 语 指导教师： 年 月 日 Broyden方法求解非线性方程组的Matlab实现 注：matlab代码来自网络，仅供学习参考。 1.把以下代码复制在一个.m文件上 function [sol, it_hist, ierr] = brsola(x,f,tol, parms) % Broyden's Method solver, globally convergent % solver for f(x) = 0, Armijo rule, one vector storage % % This code es with no guarantee or warranty of any kind. % % function [sol, it_hist, ierr] = brsola(x,f,tol,parms) % % inputs: % initial iterate = x % function = f % tol = [atol, rtol] relative/absolute % error tolerances for the nonlinear iteration % parms = [maxit, maxdim] % maxit = maxmium number of nonlinear iterations % default = 40 % maxdim = maximum number of Broyden iterations % before restart, so maxdim-1 vectors are % stored % default = 40 % % output: % sol = solution % it_hist(maxit,3) = scaled l2 norms of nonlinear residuals % for the iteration, number function evaluations, % and number of steplength reductions % ierr = 0 upon successful termination % ierr = 1 if after maxit iterations % the termination criterion is not satsified. % ierr = 2 failure in the line search. The iteration % is terminated if too many steplength reductions % are taken. % % % internal parameter: % debug = turns on/off iteration statistics display as % the iteration progresses 解线性方程组的迭代法 1.rs里查森迭代法求线性方程组Ax=b的解 function[x,n]=rs(A,b,x0,eps,M) if(nargin==3) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值elseif(nargin==4) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-A)*x0+b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 2.crs里查森参数迭代法求线性方程组Ax=b的解 function[x,n]=crs(A,b,x0,w,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1; %迭代过程 while(tol>eps) x=(I-w*A)*x0+w*b; n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 3.grs里查森迭代法求线性方程组Ax=b的解 function[x,n]=grs(A,b,x0,W,eps,M) if(nargin==4) eps=1.0e-6;%eps表示迭代精度 M=10000;%M表示迭代步数的限制值 elseif(nargin==5) M=10000; end I=eye(size(A)); n=0; x=x0; tol=1;%前后两次迭代结果误差 %迭代过程 while(tol>eps) x=(I-W*A)*x0+W*b;%迭代公式 n=n+1;%n为最终求出解时的迭代步数tol=norm(x-x0); x0=x; if(n>=M) disp('Warning:迭代次数太多，可能不收敛！'); return; end end 4.jacobi雅可比迭代法求线性方程组Ax=b的解 function[x,n]=jacobi(A,b,x0,eps,varargin) if nargin==3 eps=1.0e-6; M=200; elseif nargin<3 error return elseif nargin==5 M=varargin{1}; end D=diag(diag(A));%求A的对角矩阵 L=-tril(A,-1);%求A的下三角阵 MatLab解线性方程组一文通 当齐次线性方程AX=0,rank(A)=r 在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法. 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB 求解线性方程组 solve，linsolve 例： A=[5 0 4 2;1 -1 2 1;4 1 2 0;1 1 1 1]; %矩阵的行之间用分号隔开，元素之间用逗号或空格 B=[3;1;1;0] X=zeros(4,1);%建立一个4元列向量 X=linsolve(A,B) diff（fun，var，n）：对表达式fun中的变量var求n阶导数。 例如：F=sym（'u(x,y)*v(x,y)'）; %sym（）用来定义一个符号表达式 diff(F); %matlab区分大小写 pretty(ans) %pretty（）：用习惯书写方式显示变量；ans是答案表达式 非线性方程求解 fsolve(fun,x0,options) 为待解方程或方程组的文件名；fun其中 x0位求解方程的初始向量或矩阵； option为设置命令参数 建立文件fun.m： function y=fun(x) y=[x(1)-0.5*sin(x(1))-0.3*cos(x(2)), ... x(2) - 0.5*cos(x(1))+0.3*sin(x(2))]; >>clear;x0=[0.1,0.1];fsolve(@fun,x0,optimset('fsolve')) 注： ...为续行符 m文件必须以function为文件头，调用符为@；文件名必须与定义的函数名相同；fsolve（）主要求解复杂非线性方程和方程组，求解过程是一个逼近过程。Matlab求解线性方程组 AX=B或XA=B 在MATLAB中，求解线性方程组时，主要采用前面章节介绍的除法运算符“/”和“\”。如： X=A\B表示求矩阵方程AX＝B的解； 的解。XA=B表示矩阵方程B/A＝X． 对方程组X＝A\B，要求A和B用相同的行数，X和B有相同的列数，它的行数等于矩阵A的列数，方程X＝B/A同理。 如果矩阵A不是方阵，其维数是m×n，则有： m＝n 恰定方程，求解精确解； m>n 超定方程，寻求最小二乘解； m 中南大学MATLAB程序设计实践学长有爱奉献，下载填上信息即可上交，没有下载券的自行百度。所需m文件照本文档做即可，即新建（FILE）→脚本（NEW-Sscript）→复制本文档代码→运行（会跳出保存界面，文件名默认不要修改，保存）→结果。第一题需要把数据文本文档和m文件放在一起。全部测试无误，放心使用。本文档针对做牛顿法求非线性函数题目的同学，当然第一题都一样，所有人都可以用。←记得删掉这段话 班级： ？ 学号： 姓名： 一、《MATLAB程序设计实践》Matlab基础 表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散 空间函数值来表示取向分布函数，是三维取向分布函数的一个实例。 由于数据量非常大，不便于分析，需要借助图形来分析。请你编写一 个matlab程序画出如下的几种图形来分析其取向分布特征： （1）用Slice函数给出其整体分布特征； " ~ （2）用pcolor或contour函数分别给出(φ2＝0, 5, 10, 15, 20, 25, 30, 35 … 90)切面上f分布情况(需要用到subplot函数)； (3) 用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。 ( 备注：数据格式说明 解： （1）（ （2）将文件内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如 下： fid=fopen(''); for i=1:18 tline=fgetl(fid); end phi1=1;phi=1;phi2=1;line=0; f=zeros(19,19,19); [ while ~feof(fid) tline=fgetl(fid); data=str2num(tline); line=line+1;数据说明部分，与 作图无关此方向表示f随着 φ1从0,5,10,15, 20 …到90的变化而 变化 此方向表示f随着φ 从0,5,10,15, 20 … 到90的变化而变化 表示以下数据为φ2＝0的数据，即f（φ1，φ，0） 3.1 方程组的逆矩阵解法及其MATLAB 程序 3.1.3 线性方程组有解的判定条件及其MATLAB 程序 判定线性方程组A n m ?b X =是否有解的MATLAB 程序 function [RA,RB,n]=jiepb(A,b) B=[A b];n=length(b); RA=rank(A); RB=rank(B);zhica=RB-RA; if zhica>0, disp('请注意：因为RA~=RB ，所以此方程组无解.') return end if RA==RB if RA==n disp('请注意：因为RA=RB=n ，所以此方程组有唯一解.') else disp('请注意：因为RA=RB 遗传算法解非线性方程组的Matlab程序 程序用MATLAB语言编写。之所以选择MATLB，是因为它简单，但又功能强大。写1行MATLAB程序，相当于写10行C++程序。在编写算法阶段，最好用MATLAB语言，算法验证以后，要进入工程阶段，再把它翻译成C++语言。 本程序的算法很简单，只具有示意性，不能用于实战。 非线性方程组的实例在函数(2)nonLinearSumError1(x)中，你可以用这个实例做样子构造你自己待解的非线性方程组。 %注意：标准遗传算法的一个重要概念是，染色体是可能解的2进制顺序号，由这个序号在可能解的集合(解空间)中找到可能解 %程序的流程如下： %程序初始化，随机生成一组可能解(第一批染色体) %1: 由可能解的序号寻找解本身(关键步骤) %2：把解代入非线性方程计算误差，如果误差符合要求，停止计算 %3：选择最好解对应的最优染色体 %4：保留每次迭代产生的最好的染色体，以防最好染色体丢失 %5: 把保留的最好的染色体holdBestChromosome加入到染色体群中 %6: 为每一条染色体(即可能解的序号)定义一个概率(关键步骤) %7：按照概率筛选染色体(关键步骤) %8：染色体杂交(关键步骤) %9：变异 %10:到1 %这是遗传算法的主程序，它需要调用的函数如下。 %由染色体(可能解的2进制)顺序号找到可能解： %(1)x=chromosome_x(fatherChromosomeGroup,oneDimensionSet,solutionSum); %把解代入非线性方程组计算误差函数：(2)functionError=nonLinearSumError1(x); %判定程是否得解函数：(3)[solution,isTrue]=isSolution(x,funtionError,solutionSumError); %选择最优染色体函数： %(4)[bestChromosome,leastFunctionError]=best_worstChromosome(fatherChromosomeGroup,functionError); %误差比较函数：从两个染色体中，选出误差较小的染色体 %(5)[holdBestChromosome,holdLeastFunctionError]... % =compareBestChromosome(holdBestChromosome,holdLeastFunctionError,... % bestChromosome,leastFuntionError) %为染色体定义概率函数，好的染色体概率高，坏染色体概率低 %(6)p=chromosomeProbability(functionError); %按概率选择染色体函数： %(7)slecteChromosomeGroup=selecteChromome(fatherChromosomeGroup,p); %父代染色体杂交产生子代染色体函数 %(8)sonChrmosomeGroup=crossChromosome(slecteChromosomeGroup,2); %防止染色体超出解空间的函数 %(9)chromosomeGroup=checkSequence(chromosomeGroup,solutionSum) %变异函数 %(10)fatherChromosomeGroup=varianceCh(sonChromosomeGroup,0.8,solutionN); %通过实验有如下结果： %1。染色体应当多一些 实验1.1 用matlab 求解线性方程组 第一节 线性方程组的求解 一、齐次方程组的求解 rref （A ） %将矩阵A 化为阶梯形的最简式 null （A ） %求满足AX ＝0的解空间的一组基，即齐次线性方程组的基 础解系 【例１】 求下列齐次线性方程组的一个基础解系，并写出通解： 我们可以通过两种方法来解： 解法1： >> A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; >> rref(A) 执行后可得结果： ans= 1 -1 0 0 0 0 -1 1 0 0 0 0 由最简行阶梯型矩阵，得化简后的方程 ??? ??=+--=+--=-+-0 22004321 43214321x x x x x x x x x x x x 取x2，x4为自由未知量，扩充方程组为 即 提取自由未知量系数形成的列向量为基础解系，记 所以齐次方程组的通解为 解法2： clear A=[1 -1 1 -1;1 -1 -1 1;1 -1 -2 2]; B=null(A, 'r') % help null 看看加个‘r’是什么作用, 若去掉r ，是什么结果？ 执行后可得结果： B= 1 0 1 0 0 1 0 1 ?? ?=-=-0 04321x x x x ?????? ?====4 4432221x x x x x x x x ??? ??? ??????+????????????=????? ???????1100001142 4321x x x x x x , 00111????? ? ??????=ε, 11002????? ???????=ε2 211εεk k x += 3-7变量非线性方程组的逆Broyden解法 matlab程序 function n_broyden clear all %清内存 clc %清屏 format long i=input('请输入非线性方程组个数(3-7)i='); switch i case 3 syms x y z x0=input('请输入初值（三维列向量[x;y;z])='); eps=input('请输入允许的误差精度eps='); f1=input('请输入f1(x,y,z)='); f2=input('请输入f2(x,y,z)='); f3=input('请输入f3(x,y,z)='); F=[f1;f2;f3]; J=jacobian(F,[x,y,z]); %求jacobi矩阵 Fx0=subs(F,{x,y,z},x0); Jx0=subs(J,{x,y,z},x0); H=inv(Jx0); x1=x0-H*Fx0; k=0; while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0; q=subs(F,{x,y,z},x1)-subs(F,{x,y,z},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1; x0=x1; Fx0=subs(F,{x,y,z},x0); x1=x1-H*Fx0; k=k+1; end x1 k case 4 syms a b c d x0=input('请输入初值（列向量[a;b;c;d])='); eps=input('请输入允许的误差精度eps='); f1=input('请输入f1(a,b,c,d)='); f2=input('请输入f2(a,b,c,d)='); f3=input('请输入f3(a,b,c,d)='); f4=input('请输入f4(a,b,c,d)='); F=[f1;f2;f3;f4]; J=jacobian(F,[a,b,c,d]); %求jacobi矩阵 Fx0=subs(F,{a,b,c,d},x0); Jx0=subs(J,{a,b,c,d},x0); H=inv(Jx0); x1=x0-H*Fx0; k=0; while norm(x1-x0)>eps %用2范数来做循环条件 p=x1-x0; q=subs(F,{a,b,c,d},x1)-subs(F,{a,b,c,d},x0); H=H-((H*q-p)*p'*H)*(p'*H*q)^-1; x0=x1; Fx0=subs(F,{a,b,c,d},x0); x1=x1-H*Fx0; k=k+1; end x1 k case 5 syms a b c d e x0=input('请输入初值（列向量[a;b;c;d;e])='); eps=input('请输入允许的误差精度eps='); f1=input('请输入f1(a,b,c,d,e)='); f2=input('请输入f2(a,b,c,d,e)='); f3=input('请输入f3(a,b,c,d,e)='); f4=input('请输入f4(a,b,c,d,e)='); f5=input('请输入f5(a,b,c,d,e)='); F=[f1;f2;f3;f4;f5]; J=jacobian(F,[a,b,c,d,e]); %求jacobi矩阵 Fx0=subs(F,{a,b,c,d,e},x0); Jx0=subs(J,{a,b,c,d,e},x0); H=inv(Jx0); x1=x0-H*Fx0; 第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 x a=[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 的导函数 Matlab线性方程组求解 1. Gauss消元法： function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k); %计算行列式 end det=det*a(n,n); for k=n:-1:1 %回代求解 for j=k+1:n b(k)=b(k)-a(k,j)*x(j); end x(k)=b(k)/a(k,k); end Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 2. 列主元Gauss消去法： function x=detGauss(a,b) % Gauss列主元消去法 [n,m]=size(a); nb=length(b); det=1; %存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0; %选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n 中南大学 MATLAB程序设计实践学长有爱奉献，下载填上信息即可上交，没有下载券 的自行百度。所需m文件照本文档做即可，即新建（FILE）→脚本（NEW-Sscript）→复制本文档代码→运行（会跳出 保存界面，文件名默认不要修改，保存）→结果。第 一题需要把数据文本文档和m文件放在一起。全部测 试无误，放心使用。本文档针对做牛顿法求非线性函 数题目的同学，当然第一题都一样，所有人都可以用。 ←记得删掉这段话 班级： 学号： 姓名： 一、《MATLAB程序设计实践》Matlab基础 表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散空间函数值来表示取向分布函数，是三维取向分布函数的一个实例。由于数据量非常大，不便于分析，需要借助图形来分析。请你编写一个matlab程序画出如下的几种图形来分析其取向分布特征：（1）用Slice函数给出其整体分布特征； （2）用pcolor或contour函数分别给出(φ2＝0, 5, 10, 15, 20, 25, 30, 35 … 90)切面上f分布情况(需要用到subplot函数)； (3) 用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。 备注：数据格式说明 解： （1）将文件内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如下： fid=fopen(''); for i=1:18 tline=fgetl(fid); end phi1=1;phi=1;phi2=1;line=0; f=zeros(19,19,19); while ~feof(fid) tline=fgetl(fid); data=str2num(tline); line=line+1; if mod(line,20)==1 phi2=(data/5)+1; phi=1; 数据说明部分，与作图无关 此方向表示f 随着φ1从0,5,10,15, 20 …到90的变化而变化 此方向表示f 随着φ从0,5,10,15, 20 …到90的变化而变化 表示以下数据为φ2＝0的数据，即f （φ1，φ，0） 基于Matlab 实现牛顿迭代法解非线性方程组 已知非线性方程组如下 2211221212 10801080x x x x x x x ?-++=??+-+=?? 给定初值0(0,0)T x =，要求求解精度达到0.00001 首先建立函数F(x)，方程组编程如下，将F.m 保存到工作路径中： function f=F(x) f(1)=x(1)^2-10*x(1)+x(2)^2+8; f(2)=x(1)*x(2)^2+x(1)-10*x(2)+8; f=[f(1) f(2)]; 建立函数DF(x),用于求方程组的Jacobi 矩阵，将DF.m 保存到工作路径中： function df=DF(x) df=[2*x(1)-10,2*x(2);x(2)^2+1,2*x(1)*x(2)-10]; 编程牛顿迭代法解非线性方程组，将newton.m 保存到工作路径中： clear; clc x=[0,0]'; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',0,x(1),x(2)); N=4; for i=1:N y=df\f'; x=x-y; f=F(x); df=DF(x); fprintf('%d %.7f %.7f\n',i,x(1),x(2)); if norm(y)<0.0000001 break ; else end end 运行结果如下： 0 0.0000000 0.0000000 1 0.8000000 0.8800000 2 0.9917872 0.9917117 3 0.9999752 0.9999685 4 1.0000000 1.0000000 线性方程组求解 1.直接法 Gauss消元法： function x=DelGauss(a,b) % Gauss消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 for i=k+1:n if a(k,k)==0 return end m=a(i,k)/a(k,k); for j=k+1:n a(i,j)=a(i,j)-m*a(k,j); end b(i)=b(i)-m*b(k); end det=det*a(k,k); end det=det*a(n,n); for k=n:-1:1 %回代 for j=k+1:n b(k)=b(k)-a(k,j)*x(j); end x(k)=b(k)/a(k,k); end Example: >> A=[1.0170 -0.0092 0.0095;-0.0092 0.9903 0.0136;0.0095 0.0136 0.9898]; >> b=[1 0 1]'; >> x=DelGauss(A,b) x = 0.9739 -0.0047 1.0010 列主元Gauss消去法： function x=detGauss(a,b) % Gauss列主元消去法 [n,m]=size(a); nb=length(b); det=1;%存储行列式值 x=zeros(n,1); for k=1:n-1 amax=0;% 选主元 for i=k:n if abs(a(i,k))>amax amax=abs(a(i,k));r=i; end end if amax<1e-10 return; end if r>k %交换两行 for j=k:n z=a(k,j);a(k,j)=a(r,j);a(r,j)=z; end z=b(k);b(k)=b(r);b(r)=z;det=-det; end 本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（M üller ）法和迭代法的加速等）及其MATLAB 程序，求解非线性方程组的方法及其MATLAB 程序. 2.1 方程（组）的根及其MATLAB 命令 2.1.2 求解方程（组）的solve 命令 求方程f (x )=q (x )的根可以用MATLAB 命令: >> x=solve('方程f(x)=q(x)',’待求符号变量x ’) 求方程组f i (x 1,…,x n )=q i (x 1,…,x n ) (i =1,2,…,n )的根可以用MATLAB 命令: >>E1=sym('方程f1(x1,…,xn)=q1(x1,…,xn)'); ……………………………………………………. En=sym('方程fn(x1,…,xn)=qn(x1,…,xn)'); [x1,x2,…,xn]=solve(E1,E2,…,En, x1,…,xn) 2.1.3 求解多项式方程（组）的roots 命令 如果)(x f 为多项式，则可分别用如下命令求方程0)(=x f 的根，或求导数)('x f （见表 2-1）. 2.1.4 求解方程（组）的fsolve 命令 如果非线性方程（组）是多项式形式，求这样方程（组）的数值解可以直接调用上面已经介绍过的roots 命令.如果非线性方程（组）是含有超越函数，则无法使用roots 命令，需要调用MATLAB 系统中提供的另一个程序fsolve 来求解.当然，程序fsolve 也可以用于多项式方程（组），但是它的计算量明显比roots 命令的大. fsolve 命令使用最小二乘法（least squares method ）解非线性方程（组） (F X =)0 的数值解，其中X 和F （X ）可以是向量或矩阵.此种方法需要尝试着输入解X 的初始值（向量或矩阵）X 0，即使程序中的迭代序列收敛，也不一定收敛到(F X =)0的根（见例2.1.8）. fsolve 的调用格式： X=fsolve(F,X0) 输入函数)(x F 的M 文件名和解X 的初始值（向量或矩阵）X 0，尝试着解方程（组） 实用数值方法(Matlab) 综述报告题目：非线性方程（组）求根问题 小组成员 许多数学和物理问题归结为解函数方程f(x)=0。方程f(x)=0的解称为方程的根。对于非线性方程，在某个范围内往往不止一个根，而且根的分布情况可能很复杂，面对这种情况，通常先将考察的范围花费为若干个子段，然后判断哪些子段内有根，然后再在有根子段内找出满足精度要求的近似根。为此适当选取有根子段内某一点作为根的初始值近似，然后运用迭代方法使之足部精确化。这就是方程求根的迭代法。下面介绍书上的几种方法： 1、二分法 （1）方法概要： 假定函数f(x)在[a,b]上连续，且f(a)f(b)=0,则方程f(x)=0在[a,b]内一定有实根。取其中 将其二分，判断所求的根在的左侧还是右侧，得到一个新的有根区间 点 [],长度为[a,b]的一半。对新的有根区间继续实行上述二分手段，直至二分k次后有根区间[]长度 可见，如果二分过程无限继续下去，这些有限根区间最终必收敛于一点,该点就是所求的根。在实际计算过程中不可能完成这个无限过程，允许有一定的误差，则二分k+1次后 只要有根区间[]的长度小于，那么结果关于允许误差就能“准确”地满足方程f(x)=0。 (2)计算框图： 2、开方法 对于给定,求开方值 为此，可以运用校正技术设计从预报值生成校正值的迭代公式。自然希望校正值 能更好满足所给方程： 这是个关于校正量的近似关系式，如果从中删去二次项，即可化归为一次方程 解之有 从而关于校正值有如下开方公式 上述演绎过程表明，开方法的设计思想是逐步线性化，即将二次方程的求解画归为一次方程求解过程的重复。开方公式规定了预报值与校正值之间的一种函数关系,这里 为开方法的迭代函数。 3、Newton法 （1）方法概要 非线性方程求解 摘要：利用matlab软件编写程序，分别采用二分法、牛顿法和割线法求解非线性方程， 0 2= -x e x 的根，要求精确到三位有效数字，其中对于二分法，根据首次迭代结果，事先估计迭代次数，比较实际迭代次数与估计值是否吻合。并将求出的迭代序列用表格表示。对于牛顿法和割线法，至少取3组不同的初值，比较各自迭代次数。将每次迭代计算值求出，并列于表中。 关键词：matlab、二分法、牛顿法、割线法。 引言： 现实数学物理问题中，很多可以看成是解方程的问题，即f（x）=0的问题，但是除了极少简单方程的根可以简单解析出来。大多数能表示成解析式的，大多数不便于计算，所以就涉及到算法的问题，算法里面，具体求根时，一般先寻求根的某一个初始近似值，然后再将初始近似值逐步加工成满足精度要求为止，但是，我们知道，人为计算大大的加重了我们的工作量，所以大多用计算机编程，这里有很多可以计算的软件，例如matlab等等。 正文： 一、二分法 1 二分法原理：对于在区间[，]上连续不断且满足·＜0的函数， 通过不断地把函数的零点所在的区间一分为二，使区间的两个端点逐步逼近零点，进而得到零点近似值的方法叫做二分法。 2 二分法求根步骤：(1)确定区间，，验证·＜0，给定精确度；（2)求区间，的中点；(3)计算。若=，则就是函数的零点；若· ＜0，则令=；若·＜0，则令=。(4)判断是否达到精确度；即若 ＜，则得到零点近似值（或）；否则重复步骤2-4． 3 二分法具体内容：精度要求为5e-6，,解得实际迭代次数与估计值基本吻合，迭代如下表。n=2 c=0.000000 fc=-1.000000 n=11 c=-0.705078 fc=0.003065 n=3 c=-0.500000 fc=-0.356531 n=12 c=-0.704102 fc=0.001206 n= 4 c=-0.750000 fc=0.090133 n=13 c=-0.703613 fc=0.000277 n= 5 c=-0.625000 fc=-0.14463 6 n=14 c=-0.703369 fc=-0.00018 7 n=6 c=-0.687500 fc=-0.030175 n=15 c=-0.703491 fc=0.000045 n=7 c=-0.718750 fc=0.029240 n=16 c=-0.703430 fc=-0.000071 n= 8 c=-0.703125 fc=-0.000651 n=17 c=-0.703461 fc=-0.000013 n= 9 c=-0.710938 fc=0.014249 n=18 c=-0.703476 fc=0.000016推荐-Broyden方法求解非线性方程组的Matlab实现 精品
MATLAB代码 解线性方程组的迭代法
MatLab求解线性方程组
MATLAB解线性方程组的直接方法
Matlab求解线性方程组非线性方程组
matlab程序设计实践-牛顿法解非线性方程
线性方程组求解matlab实现
遗传算法解非线性方程组的Matlab程序
实验一用matlab求解线性方程组
3-7变量非线性方程组的逆Broyden解法matlab程序
MATLAB应用 求解非线性方程
Matlab线性方程组求解(Gauss消去法)
matlab程序设计实践-牛顿法解非线性方程
基于Matlab的牛顿迭代法解非线性方程组
线性方程组求解Matlab程序(精.选)
第二章非线性方程(组)的数值解法的matlab程序
MATLAB 非线性方程(组)求根
用matlab对非线性方程求解
相关主题
文本预览