matlab非线性方程求解

非线性方程的解法

1引言

数学物理中的许多问题归结为解函数方程的问题，即，0)(=x f (1.1)

这里，)(x f 可以是代数多项式，也可以是超越函数。若有数*x 为方程0)(=x f 的根，或称函数)(x f 的零点。

设函数)(x f 在],[b a 内连续，且0)()(

如何寻求根的初始值呢？简单述之，为了明确起见，不妨设)(x f 在区间],[b a 内有一个实的单根，且0)(,0)(>

0)(*)(≤+h x f x f k k (1.2)

那么所求的根必在),(h x x k k +内，这时可取k x 或h x k +作为根的初始近似值。这种方法通常称为“定步长搜索法”。另外，还是图解法、近似方程法和解析法。

2迭代法

2.1迭代法的一般概念

迭代法是数值计算中一类典型方法，不仅用于方程求根，而且用于方程组求解，矩阵求特征值等方面。迭代法的基本思想是一种逐次逼近的方法。首先取一个精糙的近似值，然后

用同一个递推公式，反复校正这个初值，直到满足预先给定的精度要求为止。

对于迭代法，一般需要讨论的基本问题是：迭代法的构造、迭代序列的收敛性天收敛速度以及误差估计。这里，主要看看解方程迭代式的构造。

对方程(1.1)，在区间],[b a 内，可改写成为：

)(x x ?=(2.1)

取],[0b a x ∈，用递推公式：

)(1k k x x ?=+, ,2,1,0=k (2.2)

可得到序列：

∞==0210}{,,,,k k k x x x x x (2.3)

当∞→k 时，序列∞=0}{k k x 有极限x ~，且)(x ?在x ~附近连续，则在式(2.2)两边极限，得，

)

~(~x x ?=即，x ~

为方程(2.1)的根。由于方式(1.1)和方程(2.1)等价，所以，x

x ~*=即，

*

lim x x k k =∞

→式(2.2)称为迭代式，也称为迭代公式；)(x ?可称为迭代函数。称求得的序列∞=0}{k k x 为迭代序列。

2.2程序和实例

下面是基于MATLAB 的迭代法程序，用迭代格式)(1n n x g p =+，求解方程)(x g x =，其中初始值为0p 。

**************************************************************************function[p,k,err,P]=fixpt(f1021,p0,tol,max1)%f1021是给定的迭代函数。%p0是给定的初始值。

%tol 是给定的误差界。%max1是所允许的最大迭代次数。%k 是所进行的迭代次数加1。%p 是不动点的近似值。%err 是误差。%P ={p1,p2,…,pn}

P(1)=p0;for k =2:max1P(k)=feval('f1021',P(k-1));k,err =abs(P(k)-P(k-1))p =P(k);if(err

****************************************************************************例2.1

用上述程序求方程0sin 2=-x x 的一个近似解，给定初始值5.00=x ，误差界为

510-。解：先用m 文件先定义一个名为f1021.m 的函数文件。

function y =f1021(x)

y =sin(x)/x;

建立一个主程序prog1021.m

clc clear all fixpt('f1021',0.5,10^(-5),20)

然后在MATLAB 命令窗口运行上述主程序，即：

>>prog1021

计算结果如下。

k =

2err =0.4589

k =

3err =

0.1052k =

4err =

0.0292k =

5err =

0.0078k =

6err =

0.0021k =

7err =

5.7408e-004k =

8err =

1.5525e-004k =

9err =

4.1975e-005k =

10err =

1.1350e-005k =

11err = 3.0688e-006

P =Columns 1through 6

0.50000.95890.85370.88290.87510.8772

Columns 7through 11

0.8766

0.87680.87670.87670.8767ans =0.8767

3二分法

3.1二分法原理

二分法是方程求解最直观、最简单的方法。二分法以连续函数的介值定理为基础的。由介值定理知道，若函数)(x f 区间],[b a 上连续，且0)(*)( b f a f ，即)(a f 和)(b f 负号相反，

则)(x f 在],[b a 内一定有实根。二分法的基本思想是：用对分区间的方法根据分点处函数)(x f 的符号逐步将有限区间缩小，使在足够小的区间内，方程有且仅有一根。下面简述其基本步骤。

首先记b b a a ==00,。用中点2

000b a x +=将区间],[00b a 等分成2个小区间：],[00x a 和],[00b x 。然后分析可能存在的三种情况：

如果0)(*)(0=x f a f ，则0x 是零点，也就是方程的根。

如果0)(*)(0

如果0)(*)(0

对有根的新区间施行同样的操作，于是得到一系列有空的区间：

],[],[],[],[221100k k b a b a b a b a ???? (3.1)

其中每1个区间的长度都是前一区间长度的一半，最后1个区间的长度为：

k k k a

b a b 2-=-(3.2)

如果取最后1个区间],[k k b a 的中点：

2k k k a b x -=(3.3)

作为0)(=x f 根的近似值，则有误差估计式：

1*22+-≤-≤-k k k k a b a b x x (3.4)

对于所给精度ε，若取k 使得

ε≤-+12k a b (3.5)则有，

ε≤-k x x *(3.6)

3.2程序与实例

用二分法求解方程0)(=x f 在有根区间],[b a 内的一个根，其中)(x f 在],[b a 只有一个根的情形。

*************************************************************************** function[c,err,yc]=bisect(f1031,a,b,delta)

%f1031是所要求解的函数。

%a和b分别为有根区间的左右限。

%delta是允许的误差界。

%c为所求近似解。

%yc为函数f在c上的值。

%err是c的误差估计。

ya=feval('f1031',a);

yb=feval('f',b);

if yb==0,

c=b；

return

end

if ya*yb>0,

disp('(a,b)不是有根区间');

return

end

max1=1+round((log(b-a)-log(delta))/log(2));

for k=1:max1

c=(a+b)/2;

yc=feval('f',c);

if yc==0

a=c;

b=c;

return

elseif yb*yc>0

b=c;

yb=yc;

else

a=c;

ya=yc;

end

if(b-a)

end

k;

c=(a+b)/2;

err=abs(b-a);

yc=feval('f',c);

*****************************************************************************

例3.1用上述程序求方程013=--x x 在区间]5.1,0.1[内的一个近似解，要求准确到小数点后第2位。

解：先用m 文件先定义一个名为f1031.m 的函数文件。

function y =f1031(x)

y =x^3-x -1;

建立一个主程序prog1031.m

clc clear all bisect('f1031',1.0,1.5,0.0005)

然后在MATLAB 命令窗口运行上述主程序，即：

>>prog1031

计算结果如下。

ans =

1.3247

4牛顿法

4.1牛顿法原理

从前面迭代法，我们知道，迭代函数)(x ?构造的好坏，不仅影响收敛速度，而且迭代格式有可能发散。怎样选择一个迭代函数才能保证迭代序列一定收敛呢？

构代迭代函数的一条重要途径是用近似方程来代替原方程去求根。因此，如果能将非线性方程(1.1)用线性方程去代替，那么，求近似根问题就很容易解决，而且十分方便。牛顿(Newton)法就是一种将非线性方程线化的一种方法。

设k x 是方程(1.1)的一个近似根，把如果)(x f 在k x 处作一阶Taylor 展开，即：

))((')()(k k k x x x f x f x f -+≈(4.1)

于是我们得到如下近似方程：

0))((')(=-+k k k x x x f x f (4.2)

设0)('≠k x f ，则方程(10.2.1)的解为：

()

'()k k k f x x x f x =+ (4.3)

取x ~

作为原方程(1.1)的新近似根1+k x ，即令：)(1k k k k x f x x -=+， ,2,1,0=k (4.4)

上式称为牛顿迭代格式。用牛顿迭代格式求方程的根的方法就称为牛顿迭代法，简称牛顿法。

牛顿法具有明显的几何意义。方程：

))((')(k k k x x x f x f y -+=(10.4.5)

是曲线)(x f y =上点))(,(k k x f x 处的切线方程。迭代格式(4.4)就是用切线式(4.5)的零点来代替曲线(1.1)的零点。正因为如此，牛顿法也称为切线法。

牛顿迭代法对单根至少是二阶局部收敛的，而对于重根是一阶局部收敛的。一般来说，牛顿法对初值0x 的要求较高，初值足够靠近*x 时才能保证收敛。若要保证初值在较大范围内收敛，则需对)(x f 加一些条件。如果所加的条件不满足，而导致牛顿法不收敛时，则需对牛顿法作一些改时，即可以采用下面的迭代格式：

)(')(1k k k k x f x f x x λ-=+， ,2,1,0=k (4.6)

式中，10<<λ，称为下山因子。因此，用这种方法求方程的根，也称为牛顿下山法。

牛顿法对单根收敛速度快，但每迭代一次，除需计算)(k x f 之外，还要计算)('k x f 的值。如果)(x f 比较复杂，计算)('k x f 的工作量就可能比较大。为了避免计算导数值，我们可用差商来代替导数。通常用如下几种方法：

(1)割线法。如果用

1

1)

()(----k k k k x f x f 代替)('k x f ，则得到割线法的迭代格式为：

)(111k k k k k k k x f x x x x --+---=(4.7)

(2)拟牛顿法。如果用

))

(()(1k k k k x f x f x f ---代替)('k x f ，则得到拟牛顿法的迭代格式为：

))(()()(121-+---=k k k k k k x f x f x f x f x x (4.8)

(3)Steffenson 法。如果用

)

()

())((k k k k x f x f x f x f -+代替)('k x f ，则得到拟牛顿法的迭代格式为：

)())(()(21k k k k k k x f x f x f x f x x -+-=+(4.9)

4.2程序与实例

1.牛顿法的程序

给定初值0p ，用牛顿法格式)

(')(111---=k k k k p f p f p p ， ,2,1=k ，求解非线性方程0)(=x f 。**************************************************************************

function [p1,err,k,y]=newton(f1041,df1041,p0,delta,max1)%f1041是非线性函数。%df1041是f1041的微商。%p0是初始值。%delta 是给定允许误差。%max1是迭代的最大次数。%p1是牛顿法求得的方程的近似解。%err 是p0的误差估计。%k 是迭代次数。%y =f(p1)p0,feval('f1041',p0)

for

k =1:max1p1=p0-feval('f1041',p0)/feval('df1041',p0);err =abs(p1-p0);p0=p1;p1,err,k,y =feval('f1041',p1)if (err

end ****************************************************************************例4.1用上述程序求方程0233=+-x x 的一个近似解，给定初值为1.2，误差界为610-。解：先用m 文件先定义二个名为f1041.m 和df1041.m 的函数文件。

function y =f1041(x)

y =x^3–3*x +2;

function y=df1041(x)

y=3*x^2-3;

建立一个主程序prog1041.m

clear newton('f1041','df1041',1.2,10^(-6),18)

然后在MATLAB 命令窗口运行上述主程序，即：

>>prog1041

计算结果如下。

p0=

1.2000ans =

0.1280p1=

1.1030err =

0.0970k =

1y =

0.0329p1= 1.1030

err=0.0970 k=1

y=0.0329 p1= 1.0524 err=0.0507 k=2

y=0.0084 p1= 1.0524 err=0.0507 k=2

y=0.0084 p1= 1.0264 err=0.0260 k=3

y=0.0021 p1= 1.0264 err=0.0260 k=3

y=0.0021 p1= 1.0133 err=0.0131 k=4

y= 5.2963e-004 p1= 1.0133 err=0.0131 k=4

y= 5.2963e-004 p1= 1.0066 err=0.0066 k=5

y= 1.3270e-004 p1= 1.0066 err=0.0066 k=5

y= 1.3270e-004 p1= 1.0033 err=0.0033 k=6

y= 3.3211e-005 p1= 1.0033 err=0.0033 k=6

y= 3.3211e-005 p1= 1.0017 err=0.0017 k=7

y=8.3074e-006 p1= 1.0017 err=0.0017 k=7

y=8.3074e-006 p1= 1.0008 err=8.3157e-004 k=8

y= 2.0774e-006 p1= 1.0008 err=8.3157e-004 k=8

y= 2.0774e-006 p1= 1.0004

err= 4.1596e-004 k=9

y= 5.1943e-007 p1= 1.0004 err= 4.1596e-004 k=9

y= 5.1943e-007 p1= 1.0002 err= 2.0802e-004 k=10

y= 1.2987e-007 p1= 1.0002 err= 2.0802e-004 k=10

y= 1.2987e-007 p1= 1.0001 err= 1.0402e-004 k=11

y= 3.2468e-008 p1= 1.0001 err= 1.0402e-004 k=11

y= 3.2468e-008 p1= 1.0001 err= 5.2014e-005 k=12

y=8.1170e-009 p1= 1.0001 err= 5.2014e-005 k=12

p1= 1.0000 err= 2.6008e-005 k=13

y= 2.0293e-009 p1= 1.0000 err= 2.6008e-005 k=13

y= 2.0293e-009 p1= 1.0000 err= 1.3004e-005 k=14

y= 5.0732e-010 p1= 1.0000 err= 1.3004e-005 k=14

y= 5.0732e-010 p1= 1.0000 err= 6.5020e-006 k=15

y= 1.2683e-010 p1= 1.0000 err= 6.5020e-006 k=15

y= 1.2683e-010 p1= 1.0000 err= 3.2510e-006 k=16

y= 3.1708e-011 p1= 1.0000

k=16

y= 3.1708e-011

p1= 1.0000

err= 1.6255e-006

k=17

y=7.9270e-012

p1= 1.0000

err= 1.6255e-006

k=17

y=7.9270e-012

p1= 1.0000

err=8.1277e-007

k=18

y= 1.9817e-012

ans= 1.0000

这说明，经过18次迭代得到满足精度要求的值。

2.割线法的MATLAB实现

**************************************************************************** function[p1,err,k,y]=secant(f1042,p0,p1,delta,max1)

%f1042是给定的非线性函数。

%p0,p1为初始值。

%delta为给定误差界。

%max1是迭代次数的上限。

%p1为所求得的方程的近似解。

%err为p1-p0的绝对值。

%k为所需的迭代次数。

%y＝f(p1)

p0,p1,feval('f1042',p0),feval('f1042',p1),k=0;

for k=1:max1

p2=p1-feval('f1042',p1)*(p1-p0)/(feval('f1042',p1)-feval('f1042',p0));err =abs(p2-p1);p0=p1;p1=p2;p1,err,k,y =feval('f1042',p1)；if (err

end ******************************************************************************

例4.2用上述程序求方程023=+-x x 的一个近似解，给定初值为-1.5，-1.52,误差界为610-。

解：先用m 文件先定义一个名为f1042.m 的函数文件。

function y =f1042(x)

y =x^3–x +2;

建立一个主程序prog1042.m

clc clear secant('f1042',-1.5,-1.52,10^(-6),11)

然后在MATLAB 命令窗口运行上述主程序，即：

>>prog1042

计算结果如下。

p0=

-1.5000p1=

-1.5200ans =

0.1250ans =

0.0082p1=

-1.5214err =

0.0014k =

1p1=-1.5214

err= 2.2961e-005

k=2

p1=-1.5214

err= 2.4318e-008

k=3

ans=-1.5214

这就表明，经过3次迭代得到满足精度要求的的近似解-1.5214。

推荐-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-非线性规划-作业

实例matlab-非线性规划-作业

现代设计方法-工程优化理论、方法与设计 姓名 学号 班级 研 问题 ： 某厂向用户提供发动机，合同规定，第一、二、三季度末分别交货40台、60台、80台。每季度的生产费用为 （元），其中x 是该季生产的台数。若交货后有剩余，可用于下季度交货，但需支付存储费，每台每季度c 元。已知工厂每季度最大生产能力为100台，第一季度开始时无存货，设a=50、b=0.2、c=4，问工厂应如何安排生产计划，才能既满足合同又使总费用最低。讨论a 、b 、c 变化对计划的影响，并作出合理的解释。 问题的分析和假设： 问题分析：本题是一个有约束条件的二次规划问题。决策变量是工厂每季度生产的台数，目标函数是总费用（包括生产费用和存储费）。约束条件是生产合同，生产能力的限制。在这些条件下需要如何安排生产计划，才能既满足合同又使总费用最低。 问题假设： 1、工厂最大生产能力不会发生变化； 2、合同不会发生变更； 3、第一季度开始时工厂无存货； 4、生产总量达到180台时，不在进行生产； 5、工厂生产处的发动机质量有保证，不考虑退货等因素； 6、不考虑产品运输费用是否有厂家承担等和生产无关的因素。 符号规定： x1——第一季度生产的台数； x2——第二季度生产的台数； 180-x1-x2——第三季度生产的台数； y1——第一季度总费用； y2——第二季度总费用； y3——第三季度总费用； y ——总费用（包括生产费用和存储费）。 ()2bx ax x f +=

建模： 1、第一、二、三季度末分别交货40台、60台、80台； 2、每季度的生产费用为 （元）； 3、每季度生产数量满足40 ≤x1≤100，0≤x2≤100，100≤x1+x2 ≤180； 4、要求总费用最低，这是一个目标规划模型。 目标函数: y1 2111x b x a Z ?+?= y2()4012222-?+?+?=x c x b x a Z y3()()()10018018021221213 -+?+--?+--?=x x c x x b x x a Z y x x x x x x Z Z Z Z 68644.04.04.0149201 212221321--+++=++= 40≤x1≤100 0≤x2≤100 100≤x1+x2≤180 ()2 bx ax x f +=

运用Matlab进行线性规划求解(实例)

线性规划 线性规划是处理线性目标函数和线性约束的一种较为成熟的方法，目前已经广泛应用于军事、经济、工业、农业、教育、商业和社会科学等许多方面。 8.2.1 基本数学原理 线性规划问题的标准形式是： ????? ??????≥=+++=+++=++++++=0,,,min 21221122222121112 121112211n m n mn m m n n n n n n x x x b x a x a x a b x a x a x a b x a x a x a x c x c x c z 或 ???? ?????=≥===∑∑==n j x m i b x a x c z j n j i j ij n j j j ,,2,1,0,,2,1,min 1 1 写成矩阵形式为： ?? ???≥==O X b AX CX z min 线性规划的标准形式要求使目标函数最小化，约束条件取等式，变量b 非负。不符合这几个条件的线性模型可以转化成标准形式。 MATLAB 采用投影法求解线性规划问题，该方法是单纯形法的变种。 8.2.2 有关函数介绍 在MATLAB 工具箱中，可用linprog 函数求解线性规划问题。 linprog 函数的调用格式如下： ●x=linprog(f,A,b)：求解问题minf'*x ，约束条件为A*x<=b 。 ●x=linprog(f,A,b,Aeq,beq)：求解上面的问题，但增加等式约束，即Aeq*x=beq 。若没有不等式约束，则令A=[ ],b=[ ]。 ●x=linprog(f,A,b,Aeq,beq,lb,ub)：定义设计x 的下界lb 和上界ub ，使得x 始终在该范围内。若没有等式约束，令Aeq=[ ],beq=[ ]。 ●x=linprog(f,A,b,Aeq,beq,lb,ub,x0)：设置初值为x0。该选项只适用于中型问题，默认时大型算法将忽略初值。 ●x=linprog(f,A,b,Aeq,beq,lb,ub,x0,options)：用options 指定的优化参数进行最小化。 ●[x,fval]=linprog(…)：返回解x 处的目标函数值fval 。 ●[x,lambda,exitflag]=linprog(…)：返回exitflag 值，描述函数计算的退出条件。 ●[x,lambda,exitflag,output]=linprog(…)：返回包含优化信息的输出参数output 。 ●[x,fval,exitflag,output,lambda]=linprog(…)：将解x 处的拉格朗日乘子返回到lambda 参数中。

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

MATLAB解线性方程组的直接方法

在这章中我们要学习线性方程组的直接法，特别是适合用数学软件在计算机上求解的方法. 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> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7]; b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果为 请注意：因为RA=RB=n ，所以此方程组有唯一解. RA = 4,RB =4,n =4 在MATLAB 工作窗口输入 >>X=A\b, 运行后输出结果为 X =(0 0 0 0)’. (2) 在MATLAB 工作窗口输入程序 >> A=[3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3];b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b)

matlab程序设计实践-牛顿法解非线性方程

中南大学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求解线性方程组非线性方程组

求解线性方程组 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非线性规划

一般非线性规划 标准型为： min F(X) s.t AX<=b b e q X A e q =? G(X)0≤ Ceq(X)=0 VLB ≤X ≤VUB 其中X 为n 维变元向量，G(X)与Ceq(X)均为非线性函数组成的向量，其它变量的含义与线性规划、二次规划中相同.用Matlab 求解上述问题，基本步骤分三步： 1. 首先建立M 文件fun.m,定义目标函数F （X ）: function f=fun(X); f=F(X); 2. 若约束条件中有非线性约束:G(X)0≤或Ceq(X)=0,则建立M 文件nonlcon.m 定义函数G(X)与Ceq(X): function [G,Ceq]=nonlcon(X) G=... Ceq=... 3. 建立主程序.非线性规划求解的函数是fmincon,命令的基本格式如下： (1) x=fmincon (‘fun’,X0,A,b) (2) x=fmincon (‘fun’,X0,A,b,Aeq,beq) (3) x=fmincon (‘fun’,X0,A,b, Aeq,beq,VLB,VUB) (4) x=fmincon (‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’) (5)x=fmincon (‘fun’,X0,A,b,Aeq,beq,VLB,VUB,’nonlcon’,options) (6) [x,fval]= fmincon(...) (7) [x,fval,exitflag]= fmincon(...) (8)[x,fval,exitflag,output]= fmincon(...) 注意： [1] fmincon 函数提供了大型优化算法和中型优化算法。默认时，若在

线性方程组求解matlab实现

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> A=[2 3 -1 5;3 1 2 -7;4 1 -3 6;1 -2 4 -7]; b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果为 请注意：因为RA=RB=n ，所以此方程组有唯一解. RA = 4,RB =4,n =4 在MATLAB 工作窗口输入 >>X=A\b, 运行后输出结果为 X =(0 0 0 0)’. (2) 在MATLAB 工作窗口输入程序 >> A=[3 4 -5 7;2 -3 3 -2;4 11 -13 16;7 -2 1 3];b=[ 0; 0; 0; 0]; [RA,RB,n]=jiepb(A,b) 运行后输出结果 请注意：因为RA=RB> A=[4 2 -1;3 -1 2;11 3 0]; b=[2;10;8]; [RA,RB,n]=jiepb(A,B) 运行后输出结果 请注意：因为RA~=RB ，所以此方程组无解. RA =2,RB =3,n =3 (4)在MATLAB 工作窗口输入程序

遗传算法解非线性方程组的Matlab程序

遗传算法解非线性方程组的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。染色体应当多一些

遗传算法解决非线性规划问题的Matlab程序

通常，非线性整数规划是一个具有指数复杂度的NP问题，如果约束较为复杂，Matlab优 化工具箱和一些优化软件比如lingo等，常常无法应用，即使能应用也不能给出一个较为令 人满意的解。这时就需要针对问题设计专门的优化算法。下面举一个遗传算法应用于非线性整数规划的编程实例，供大家参考! 模型的形式和适应度函数定义如下: nun ￡ =迟叼匸［（1_冏）督 i-1 /-I J=K乙员-??严丿=12 M…严 ▼ 0 或1* 适应度函数为3 Fi tn叱O）=》〔?巾1口｛＞?（卡（￡）一/；0?门））转幷亠 Z j'-i 50 4 S0 其中比=2、即士￡ = ￡ =瓦％■，口（1-务），马；j^ = s = ■ x v' y- to.8,02）., /-I i-L i-1 E 这是一个具有200个01决策变量的多目标非线性整数规划，编写优化的目标函数如下，其 中将多目标转化为单目标采用简单的加权处理。 fun ctio n Fit ness=FITNESS(x,FARM,e,q,w) %%适应度函数 %输入参数列表 % x 决策变量构成的 4X50的0-1矩阵 % FARM 细胞结构存储的当前种群，它包含了个体x

% e 4 X50的系数矩阵 % q 4 X50的系数矩阵 % w 1 X50的系数矩阵 %% gamma=0.98; N=length(FARM);% 种群规模 F1=zeros(1,N); F2=zeros(1,N); for i=1:N xx=FARM{i}; ppp=(1-xx)+(1-q).*xx; F1(i)=sum(w.*prod(ppp)); F2(i)=sum(sum(e.*xx)); end ppp=(1-x)+(1-q).*x; f1=sum(w.*prod(ppp)); f2=sum(sum(e.*x)); Fitness=gamma*sum(min([sign(f1-F1);zeros(1,N)]))+(1-gamma)*sum(min([sign(f2- F2);zeros(1,N)])); 针对问题设计的遗传算法如下，其中对模型约束的处理是重点考虑的地方 function [Xp,LC1,LC2,LC3,LC4]=MYGA(M,N,Pm) %% 求解 01 整数规划的遗传算法 %% 输入参数列表

实验一用matlab求解线性方程组

实验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应用 求解非线性方程

第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线性方程组求解(Gauss消去法)

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程序设计实践-牛顿法解非线性方程

中南大学 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的牛顿迭代法解非线性方程组

基于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

完整word版matlab解非线性规划例题

关于非线性规划问题 背景： 线性规划问题，即目标函数和约束条件都是线性函数的规划问题，但在实际工作中，还常常会遇到另一类更一般的规划问题，即目标函数和约束条件中至少有一个是非线性函数问题，即非线性规划问题。 求解方法：Matlab软件 问题： 某厂向用户提供发动机，合同规定，第一、二、三季度末分别f(x)?ax?bx^2（元）交货50台、70台、90台。，每季度的生产费用为x是该季度生产的台数，若交货有剩余可用于下季度交货，但需其中c元。已知工厂每季度最大生产能力为100支付存储费，每季度每台bca=4,问工厂如何安排每台，第一季度开始时无存货，设=0.2,=50,月生产计划，才能既满足合同又使总费用最低（包括生产费用和库存费用）。 问题分析与假设： F(x)。目标函数是总费用，记为约束条件是生产合同和生产能力的限制。 x1x2台，则第三季度生产台，第二季度生产设第一季度生产(210?x1?x2)台。则： 120?x1?x2?210 50?x1?1000?x2?100 bca=4,

=0.2,=50,由． T1?50x1?0.2x1^2，第一季度生产费用 k1?4(x1?50)，剩余品存储到下一季度的费用 T2?50x2?0.2x2^2同理可得： k2?4(x1?x2?120) T3?50(210?x1?x2)?0.2(210?x1?x2)^2 建模 总费用 F(x)?T1?T2?T3?k1?k2?10300?0.2(x1^2?x2^2)?0.2(210?x1?x2)^2?4(2x1?x2?120)先建立M-文件: a=50;b=0.2;c=4; H=diag(2*b*ones(1,3));C=[a+2*c,a+c,a]; A1=[-1,0,0;-1,-1,0];b1=[-50,-120]'; A2=[1 1 1];b2=210; v1=[0 0 0]';v2=[100 100 100]'; [x,faval,exitflag,output,lambada]=quadprog(H,C,A1,b1,A2,b2,v1,v2,[]) X2=x'*H*x/2+C*x-140*c 再建立主程序: a=50;b=0.2;c=4; H=diag(2*b*ones(1,3));C=[a+2*c,a+c,a]; A1=[-1,0,0;-1,-1,0];b1=[-50,-100]'; A2=[1 1 1];b2=210; v1=[0 0 0]';v2=[100 100 100]'; [x,faval,exitflag,output,lambada]=quadprog(H,C,A1,b1,A2,b2,v1,v2,[]) X2=x'*H*x/2+C*x-140*c 运算结果： x =

MATLAB非线性优化fmincon

active-set and sqp algorithms 不接受用户提供的海塞矩阵，对拉格朗日的海塞矩阵提供一个拟牛顿的近似值； 目标函数估值次数与迭代次数？ 优化成功或失败 一、求解失败 1、在到达迭代次数阈值或目标函数估值次数阈值时，求解器没有最小化目标到要求的精度，此时求解器停止。接下来，可以尝试以下方法： （1）设置‘Display’为‘iter’，查看每步的迭代信息，这些信息包括：目标函数（Fval or f(x) or Resnorm）是否是下降的；检查约束越界（Max constraint）是否是递减趋向于0；查看一阶优化是否是递减趋向于0；查看置信域半径(Trust-region radius)是否下降趋向于一个小的值。若其中至少一种情况为是，就表示结果是不断改善的。如果结果是不断改善的，可以采取下边的措施：设置MaxIter、MaxFunEvals比默认值大的值，默认值可以在优化工具箱或求解器的函数参考页的优化表中查看；从最后计算出的点开始重新求解。如果结果没有改善，尝试以下其他的方法。（2）放松精度 如果TolX或TolFun太小，当求解器达到一个最小值时可能也不会识别到，这就会导致无限次徒劳的迭代。DiffMaxChange和DiffMinChange选项能影响求解器的改善，它们控制求导估计中有限差分的步长。 （3）从不同的初始点重新开始求解

（4）检查目标函数和约束函数的定义 举个例子，可以检查目标函数和非线性约束函数在某些特定点处返回正确的值。不可行的点不一定导致函数的错误。 （5）对问题进行中心化和标准化 当每个坐标轴对目标函数和约束函数有相同的影响时，求解器更能可靠的运行，对每个坐标轴方向乘以合适的量使得每个坐标轴的影响相同，在特定的坐标轴上加上合适的值使得它们长度一致。 （6）提供解析的梯度和雅可比矩阵 如果用户不提供解析的梯度或雅可比矩阵，求解器会用有限差分来估计这些值，因此提供这些导数可以减少运算时间，提高计算准确度。 对于约束问题，提供梯度还有另一个好处----求解器到达一个点x 时能满足该点是可行的，但有限差分在x点周围可能会导致不可行的点，在这种情况下，求解器可能会失败或突然中断。（7）提供海塞矩阵 当提供海塞矩阵时，求解器能运行的更可靠，而且运行的次数比较少。 2、无可行点 在TolCon约束精度内，求解器不能找到一个满足所有约束条件的点，此时，可以尝试以下方法： （1）检查线性约束