matlab处理非线性误差估计
用matlab拟和模型参数和计算参数误差
Matlab用以建立数学模型是一个很好的工具。对模型函数的评价，一个很重要的方法就是最小二乘（Least squares）由least mean squares这个方法得到。假如有点集P(X, Y)，每一个点 P(i) 由X(i), Y(i) , i = 1 ~ m组成；模型 Y_fit = F( A, X ), Y_fit(i) = F(A, X(i) ); 其中A= A(1) A(2) … A(n)是模型的n个参数。least mean squares = (1/m) * sum ((Y(i) - Y_fit(i) ).^2) (i = 1 ~ m)。
一个好的模型，least mean squares就小；而另一方面，如何得到模型参数A，使得least mean squares有最小值，就是所谓的，最小二乘拟合（least squares curve fitting）了。
简介：
模型有线性和非线性之分。对于线性模型，求参数，其实就是求一步矩阵的逆（稍候我们可以看到）。而非线性模型，往往不能一步就得到结果，所以就需要多步逼近。就这样，在众多的多步逼近的方法中，最快收敛于最佳参数值的方法就比较垂青。这中间，最强的当然就是Newton 法：
A: n+1 = A: n + (Hessen ( L ))^-1 * grad(L)
这里Hessen ( L )是被拟合的模型函数的least mean squares方法的Hessen
矩阵。grad(L)是她的梯度矩阵。参数矩阵A的当前值是A:n和下一步值A: n+1。
这个方法包含了一个求hessen矩阵的逆的运算。其实，这个方法难的不是这个逆，而是如何得到Hessen矩阵和梯度矩阵。梯度矩阵还好说，就是least mean squares方法的对各个参数的一介偏导数。而Hessen矩阵包含了一介偏导数的组合（主要是相乘），和二介偏导数。当然，许多模型的二介偏导数相对于一介偏导数的组合是一个比较小的量，特别是线性模型，就没有二介偏导（所以，线性模型可以直接求出参数）。于是，新的方法就利用这个特点，将逼近限制在一介偏导数构成的伪Hessen矩阵上。这就诞生了两个比较著名的方法
Gauss-Newton 法和Levenberg-Marquardt法。
Gauss-Newton 法直接用Jacobian 行列式代替 Hessen矩阵，用least squares 值代替梯度（注意，不是least mean squares，因为当用Jacobian 行列式代替Hessen矩阵时，中间有一个自由度的差别）这里的拟合就变成了
A: n+1 = A: n + (Jacobian ( L ))^-1 * L （对L的定义会在下文中给出）
因为越是接近最佳值（或者临界值），Jacobian ( L )就越是畸形，所以在实际的计算机运算中，求逆这一步都用所谓的帽子运算符 (假如 J= Jacobian
( L ) )；
( Jacobian ( L ) )^-1 ---> ( ( Trans(J) * J )^-1) * Trans(J)
这里Trans()是转置运算。
Levenberg-Marquardt法比Gauss-Newton好，可以说，在很多时候，和Newton 法是一样的（当参数互相独立的时候，而这个条件是好的数学模型基本的条件，虽然有时候不能100%保证，但是，基本上参数互相之间相关性不大）。与Newton 法相比，Levenberg-Marquardt法因为没有纯的二介偏导运算，速度有时候还比Newton 法要快。而且，Levenberg-Marquardt法引入了一个步长，使得该方法在计算机运算上，更容易控制。
实现：
下面我们简单说一下这个Levenberg-Marquardt法，和拟合过程中的参数误差估计：
假设 L = sum ((Y - Y_fit ).^2) (Y= [Y (1), Y(2)… Y(m)]是测量值，
Y_fit = F( A, X )=[ F( A(1), A(2)…A(n), X(1)), F( A(1), A(2)…A(n), X(2))…F( A(1), A(2)…A(n), X(m))] ), 是模型值。其中的X(1)~(m)就是对应于Y的测量值，A(1)~(n)是模型的参数)
使L最小的参数A就是我们要求的参数。当L线性时，L的最小值，就是要求满足L对于各个参数的偏导数=0的点（相当于解一个多元一次方程组）。当L非线性时，上述的条件依然有效，就是不能解。于是就有了逼近法。考虑到非线性模型的二介偏导数不全为零，如果给定一个接近最佳参数值的估计参数
A:current，最佳参数（使得L最小的参数）A:min可由Newton 法给出：
A:min = A:current + Hessen (L)^-1 * grad(L)。
grad(L)是一个一维向量，第j个向量g(j)是L对参数A(j)的一介偏导数，所以，有n个元素。
假设F( A, X )对A(j)的一介偏导数为D_Y_fit(j)，这样：
g(j)= - 2 * sum[(Y - Y_fit)* D_Y_fit(j)]
Hessen (L)是一个二维的矩阵，第i行第j列的值
H(i, j)=2*sum[D_Y_fit(i)* D_Y_fit(j)- (Y - Y_fit)* DD_Y_fit(j, i)]
这里DD_Y_fit(j, i)是模型函数F( A, X )先对A(j)参数求一介偏导再对A(i)参数求第二介偏导。
在Levenberg-Marquardt法中，首先考虑到(Y - Y_fit)* DD_Y_fit(j, i)这一项，在接近最佳值的时候，(Y - Y_fit)很小；其次，假如参数的互相之间相关性很小，这就使得DD_Y_fit(j, i)在i, j不相等的时候值很小甚至等于0（参数之间相互正交）。再次，考虑到2和-2是常数，可以互相抵消。于是，LM法用H’(i, j)（几何值相当于Hessen矩阵的一半）来代替Hessen矩阵，具体的计算方法如下：
H’(i, j)= sum[D_Y_fit(i)* D_Y_fit(j)] (i不等于j 时)
H’(i, i)= (1 + r) * sum[D_Y_fit(i)* D_Y_fit(i)]
假设D_A = ( A:min - A:current)是最佳参数和当前估计参数的差，是一个一
维向量。
LM法就是
1) 假设一个r比如可以从1开始，解方程D_A * H’ = (-1/2)grad(L)
细心的同学已经发现了，这是一个线性方程，意味着在matlab里面一步就可以解决：
D_A = H’ \ ((-1/2)grad(L))
没有错，但是！（见后面）
2) 比较L(A:current + D_A ) 和L(A:current )，这个时候假如
a) L(A:current + D_A ) < L(A:current )，就保留A:current + D_A，并减小r（比如r=r/10）
b) 否则就增加 r（比如r=r*10）
3) 重复解方程D_A * H’ = (-1/2)grad(L)
这个方法可以有两个结束指标a) r很大时，说明方程到了平坦区。b) 每一次解方程以后，会比较一次L，假如L还在减小，但是减小的幅度很小时，可以结束。实际应用中，没有用r很小来结束，因为L变化的减小比r快。当然，r很小也是一个结束条件，只是，r还没有很小时，L的变化已经很小了。说明拟合已经十分接近终点了。
有一点要注意，就是接近终点时计算Hessen矩阵（Newton 法）或者H’
（ Levenberg-Marquardt法）的inv会变得十分不稳定，就算matlab说用QR decomposition，也可能出现错误。如果你用r=1e+20计算时，基本上matlab backslash 是得不到什么结果的。建议用函数inv()来计算。或者用经过“最大绝对值优先”的Gauss-Jordan法。他们都比较稳定。
关于拟合的误差：
最小二乘的拟合误差可以用1/2Hessen矩阵来估计，Levenberg-Marquardt法中可以用r=0的时候的H’的inv来计算。具体的：
A的无偏标准误差（覆盖68%的置信区间）= +/- sqrt( L/(自由度)* diag(H’^-1)) 这里，L就是用最后得到的拟合参数计算的方差和，自由度的定义是
自由度 = 拟合样本数(x-y 的个数，m) - 拟合参数个数 (n)
总结：
matlab计算非线性最小二乘的方法大致就是这样一个过程：
1) 确定拟合函数F( A, X )，写出函数对各个参数的一介偏导数的函数
D_Y_fit(1~n)（如果手算有困难，可以用Mathematica之类的支持符号运算的软件帮忙，据说也可以用matlab的符号运算工具箱，不过自己没有试过），如果偏导数比较简单，还可以再写出二介偏导数DD_Y_fit(1~n, 1~n)（共n*n个，
用于Newton 法）。
2) 用一组估计值作为A:current，计算L = sum ((Y - F( A, X ) ).^2) (least square)
3) 根据1)得到的函数，用上面的定义，用估计值A:current构建(-1/2)grad(L) 矩阵G。第j个元素g'(j)可以由下面的方法得到：
g'(j)= sum[(Y - Y_fit)* D_Y_fit(j)]
和H’矩阵（Levenberg-Marquardt法）：
H’(i, j)= sum[D_Y_fit(i)* D_Y_fit(j)] (i不等于j 时)
H’(i, i)= (1 + r) * sum[D_Y_fit(i)* D_Y_fit(i)]
或者1/2Hessen矩阵（Newton 法）
H(i, j)= sum[D_Y_fit(i)* D_Y_fit(j)- (Y - Y_fit)* DD_Y_fit(j, i)]
因为A有了，X Y有了，所以上面的矩阵都是数字矩阵。对于LM法，r可以从1开始。
4)解方程
D_A= H’ \ G (LM法) 或者 D_A= H \ G (Newton 法)
用A:next = A:current + D_A 重新计算2)
注意，计算2)之前，可以对A:next作一个范围检查，对于越过边界条件的值可以修正。方法：i ) 保留原来的值；ii )用边界条件；iii)对于第j个不合条件的值，取一个小的D_A(j)（比如A:next(j) = A:current(i) + D_A(i)/10 , 重复直到满足条件）。
5) 如果是LM法，根据上面文章，改变r值，重复3)
a) L(A:next) < L(A:current )，就保留A:next (A:current = A:next)，并减小r（比如r=r/10）
b) 否则就增加 r（比如r=r*10）放弃A:next
如果是Newton 法
a) L(A:next) < L(A:current )，就保留A:next (A:current = A:next)，
b) 否则，就结束！说明A:current已经是最佳值了。
通常，Newton 法结束的时候，D_A会等于0。
6) 在重复3)之前，可以用文章提到的结束条件判断一下。
7) 完全结束了可以根据需要，求一个误差。不要的话，就结束好了。
用matlab大致就是这样。
其实写这篇文章的时候，自己正好在做关于非线性拟合的模型。发现matlab的fit不能给出参数的误差，就用以前的c++的方法改了一下，放在matlab里面。顺便总结了上面的文字，当然，只是自己的理解，如果有什么不完整的，还希望大家讨论。
下次我会帖c++的非线性拟和的方法，主要是Levenberg-Marquardt方法。大家如果有比较棘手的函数需要拟和，不妨帖一个，看看是否可以作为例子，用c++解决。
有任何问题，也欢迎大家回帖或者给自己消息！
freecat 2007-6-22 23:16
想问一下
确定拟合函数F( A, X )这一步
如果函数非线性而且非常复杂，而且不能用一个函数来描述，也就是说F未知这种情况，用什么方法求极值
[[i] 本帖最后由 freecat 于 2007-6-22 23:42 编辑 [/i]]
recbio 2007-6-23 22:13
这个问题非常好。其实自己也不是很清楚这个中间的一般方法。
通常的，如果没有规律，且非线性。如果是拟合某一个连续的区间，多项式拟合总是可以的(类似原函数的泰勒展开)。就是参数的意义不十分明确了。
如果是不连续的，那么就只有自己建立拟合规则了。
实在不行，还可以用"枚举"，建立评价函数，利用计算机的速度来解决。如果实在太多了，就用Monte Carlo法，然后看运气了。
一般来说，不能100%的解决，还是可以在一个比较高的概率范围内解题的。
当然，具体的方案还需视实际情况而定。很高兴进一步讨论。
狐狸老大 2007-6-23 23:39
建立评价函数，后用启发式的迭代优化方法求极值，是一个行之有效的方法，虽然这种方法的收敛性质未必一定稳定，但确实可以在一个比较高的概率范围内解题的。
freecat 2007-6-23 23:58
[quote]原帖由 [i]recbio[/i] 于 2007-6-23 22:13 发表
[url=http://www.dolc.de/forum/redirect.php?goto=findpost&pid=8372751& ptid=527435][img]http://www.dolc.de/forum/images/common/back.gif[/img ][/url]
这个问题非常好。其实自己也不是很清楚这个中间的一般方法。
通常的，如果没有规律，且非线性。如果是拟合某一个连续的区间，多项式拟合总是可以的(类似原函数的泰勒展开)。就是参数的意义不十分明确了。
... [/quote]
我自己的论文也是卡在了这里，对于一个所谓评价函数的求解，倒是可以用所谓启发式迭代求近似最优解，但是在函数和参数的一一对应关系不明确之前，用传统的梯度迭代法却遇到了困难。
推荐-Broyden方法求解非线性方程组的Matlab实现 精品
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
MATLAB代码 解线性方程组的迭代法
解线性方程组的迭代法 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求解线性方程组
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） 《MATLAB语言》课成论文 利用MATLAB求线性方程组 姓名：郭亚兰 学号：12010245331 专业：通信工程 班级：2010级通信工程一班 指导老师：汤全武 学院：物电学院 完成日期：2011年12月17日 利用MATLAB求解线性方程组 （郭亚兰 12010245331 2010 级通信一班） 【摘要】在高等数学及线性代数中涉及许多的数值问题，未知数的求解，微积分，不定积分，线性方程组的求解等对其手工求解都是比较复杂，而MATLAB语言正是处理线性方程组的求解的很好工具。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是数学的一个分支，它的研究对象是向量，向量空间（或称线性空间），线性变换和有限维的线性方程组。因而，线性代数被广泛地应用于抽象代数和泛函分析中；由于科学研究中的非线性模型通常可以被近似为线性模型，使得线性代数被广泛地应用于自然科学和社会科学中。线性代数是讨论矩阵理论、与矩阵结合的有限维向量空间及其线性变换理论的一门学科。 【关键字】线性代数MATLAB语言秩矩阵解 一、基本概念 1、N级行列式A：A等于所有取自不同性不同列的n个元素的积的代数和。 2、矩阵B：矩阵的概念是很直观的，可以说是一张表。 3、线性无关：一向量组（a1,a2,…,an）不线性相关，既没有不全为零的数 k1,k2,………kn使得：k1*a1+k2*a2+………+kn*an=0 4、秩：向量组的极在线性无关组所含向量的个数成为这个向量组的秩。 5、矩阵B的秩：行秩，指矩阵的行向量组的秩；列秩类似。记：R(B) 遗传算法解非线性方程组的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。染色体应当多一些 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 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 的导函数 实验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 += 中南大学 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线性方程组求解 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 实现牛顿迭代法解非线性方程组 已知非线性方程组如下 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 本章主要介绍方程根的有关概念，求方程根的步骤，确定根的初始近似值的方法（作图法，逐步搜索法等），求根的方法（二分法，迭代法，牛顿法，割线法，米勒（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）方法概要 线性方程组求解 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 非线性方程求解 摘要：利用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.000016MATLAB解线性方程组的直接方法
Matlab求解线性方程组非线性方程组
matlab程序设计实践-牛顿法解非线性方程
利用MATLAB求线性方程组
遗传算法解非线性方程组的Matlab程序
线性方程组求解matlab实现
3-7变量非线性方程组的逆Broyden解法matlab程序
MATLAB应用 求解非线性方程
实验一用matlab求解线性方程组
matlab程序设计实践-牛顿法解非线性方程
Matlab线性方程组求解(Gauss消去法)
基于Matlab的牛顿迭代法解非线性方程组
第二章非线性方程(组)的数值解法的matlab程序
MATLAB 非线性方程(组)求根
线性方程组求解Matlab程序(精.选)
用matlab对非线性方程求解
- 非线性方程组求解及matlab实现分解
- MATLAB求解非线性方程
- 非线性方程组求解及matlab实现
- 遗传算法解非线性方程组的Matlab程序
- 用matlab求解非线性方程组的几种方法之程序.
- 用Matlab求解非线性方程组-甘才军的Matlab学习笔记
- 非线性方程组求解及matlab实现(精选)
- (整理)实验2 利用matlab解(非)线性、微分方程(组).
- MATLAB 非线性方程(组)求根
- matlab非线性方程的解法(含牛拉解法)
- 实验5 非线性方程求根及其MATLAB实现资料
- MATLAB应用 求解非线性方程
- 基于Matlab的牛顿迭代法解非线性方程组
- 非线性方程组求解及matlab实现讲解
- 经典Newton-Raphson牛顿法求解非线性方程组matlab源程序
- 用matlab求解非线性方程组的几种方法之程序
- Matlab求解线性方程组非线性方程组
- 伪弧长法求解非线性方程matlab程序
- Matlab求解线性方程组、非线性方程
- matlab程序设计实践 牛顿法解非线性方程