Matlab求解线性方程组非线性方程组.docx
- 文档编号:11032189
- 上传时间:2023-02-24
- 格式:DOCX
- 页数:19
- 大小:65.62KB
Matlab求解线性方程组非线性方程组.docx
《Matlab求解线性方程组非线性方程组.docx》由会员分享,可在线阅读,更多相关《Matlab求解线性方程组非线性方程组.docx(19页珍藏版)》请在冰豆网上搜索。
Matlab求解线性方程组非线性方程组
Matlab求解线性方程组、非线性方程组
姓名:
罗宝晶学号:
1012208015专业:
材料学院高分子系
第一部分数值计算
Matlab求解线性方程组AX=B或XA=B
在MATLAB中,求解线性方程组时,主要采用除法运算符“/”和“\”。
如:
X=A\B表示求矩阵方程AX=B的解;
X=B/A表示矩阵方程XA=B的解。
对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。
如果矩阵A不是方阵,其维数是m×n,则有:
m=n恰定方程,求解精确解;
m>n超定方程,寻求最小二乘解;
m 针对不同的情况,MATLAB将采用不同的算法来求解。 一.恰定方程组 恰定方程组由n个未知数的n个方程构成,方程有唯一的一组解,其一般形式可用矩阵,向量写成如下形式: Ax=b其中A是方阵,b是一个列向量; 在线性代数中,最常用的方程组解法有: (1)利用Cramer公式来求解法; (2)利用矩阵求逆解法,即x=A-1b; (3)利用Gaussian消去法; (4)利用Lu法求解。 一般来说,对维数不高,条件数不大的矩阵,上面四种解法所得的结果差别不大。 前三种解法的真正意义是在其理论上,而不是实际的数值计算。 MATLAB中,出于对算法稳定性的考虑,行列式及逆的计算大都在Lu分解的基础上进行。 在MATLAB中,求解这类方程组的命令十分简单,直接采用表达式: x=A\b。 在MATLAB的指令解释器在确认变量A非奇异后,就对它进行Lu分解,并最终给出解x;若矩阵A的条件数很大,MATLAB会提醒用户注意所得解的可靠性。 如果矩阵A是奇异的,则Ax=b的解不存在,或者存在但不唯一;如果矩阵A接近奇异时,MATLAB将给出警告信息;如果发现A是奇异的,则计算结果为inf,并且给出警告信息;如果矩阵A是病态矩阵,也会给出警告信息。 此外还要注意: 在求解方程时,尽量不要用inv(A)*b命令,而应采用A\b的解法。 因为后者的计算速度比前者快、精度高,尤其当矩阵A的维数比较大时。 另外,除法命令的适用行较强,对于非方阵A,也能给出最小二乘解。 二.超定方程组 对于方程组Ax=b,A为n×m矩阵,如果A列满秩,且n>m。 则方程组没有精确解,此时称方程组为超定方程组。 线性超定方程组经常遇到的问题是数据的曲线拟合。 对于超定方程,在MATLAB中,利用左除命令(x=A\b)来寻求它的最小二乘解;还可以用广义逆来求,即x=pinv(A),所得的解不一定满足Ax=b,x只是最小二乘意义上的解。 左除的方法是建立在奇异值分解基础之上,由此获得的解最可靠;广义逆法是建立在对原超定方程直接进行householder变换的基础上,其算法可靠性稍逊与奇异值求解,但速度较快; 例子: 求解超定方程组 A=[2-13;31-5;4-11;13-13] A= 2-13 31-5 4-11 13-13 b=[303-6]’; rank(A) ans= 3 x1=A\b x1= 1.0000 2.0000 1.0000 x2=pinv(A)*b x2= 1.0000 2.0000 1.0000 A*x1-b ans= 1.0e-014 -0.0888 -0.0888 -0.1332 0 可见x1并不是方程Ax=b的精确解,用x2=pinv(A)*b所得的解与x1相同。 二.欠定方程组 欠定方程组未知量个数多于方程个数,但理论上有无穷个解。 MATLAB将寻求一个基本解,其中最多只能有m个非零元素。 特解由列主元qr分解求得。 例子: 解欠定方程组 A=[1-211;1-21-1;1-215] A= 1-211 1-21-1 1-21-1 1-215 b=[1-15]’ x1=A\b Warning: Rankdeficient,rank=2tol=4.6151e-015 x1= 0 -0.0000 0 1.0000 x2=pinv(A)*b x2= 0 -0.0000 0.0000 1.0000 四.方程组的非负最小二乘解 在某些条件下,所求的线性方程组的解出现负数是没有意义的。 虽然方程组可以得到精确解,但却不能取负值解。 在这种情况下,其非负最小二乘解比方程的精确解更有意义。 在MATLAB中,求非负最小二乘解常用函数nnls,其调用格式为: (1)X=nnls(A,b)返回方程Ax=b的最小二乘解,方程的求解过程被限制在x的条件下; (2)X=nnls(A,b,TOL)指定误差TOL来求解,TOL的默认值为TOL=max(size(A))*norm(A,1)*eps,矩阵的-1范数越大,求解的误差越大; (3)[X,W]=nnls(A,b)当x(i)=0时,w(i)<0;当下x(i)>0时,w(i)0,同时返回一个双向量w。 例子: 求方程组的非负最小二乘解 A=[3.4336-0.52380.6710 -0.52383.2833-0.7302 0.6710-0.73024.0261]; b=[-1.0001.50002.5000]; [X,W]=nnls(A,b) X= 0 0.6563 0.6998 W= -3.6820 -0.0000 -0.0000 x1=A\b x1= -0.3569 0.5744 0.7846 A*X-b ans= 1.1258 0.1437 -0.1616 A*x1-b ans= 1.0e-0.15 -0.2220 0.4441 0 下面再举几个用MATLAB解方程的例子 1、对于多项式p(x)=x3-6x2-72x-27,求多项式p(x)=0的根,可用多项式求根函数roots(p),其中p为多项式系数向量,即 >>p=[1,-6,-72,-27] p= 1.00-6.00-72.00-27.00 p是多项式的MATLAB描述方法,我们可用poly2str(p,'x')函数,来显示多项式的形式: >>px=poly2str(p,'x') px=x^3-6x^2-72x–27 多项式的根解法如下: >>formatrat%以有理数显示 >>r=roots(p) r= 2170/179 -648/113 -769/1980 2、在MATLAB中,求解用符号表达式表示的代数方程可由函数solve实现,其调用格式为: solve(s,v): 求解符号表达式s的代数方程,求解变量为v。 例如,求方程(x+2)x=2的解,解法如下: >>x=solve('(x+2)^x=2','x') x= .69829942170241042826920133106081 得到符号解,具有缺省精度。 如果需要指定精度的解,则: >>x=vpa(x,3) x= .698 3、使用fzero或fsolve函数,可以求解指定位置(如x0)的一个根,格式为: x=fzero(fun,x0)或x=fsolve(fun,x0)。 例如,求方程0.8x+atan(x)-π=0在x0=2附近一个根,解法如下: >>fu=@(x)0.8*x+atan(x)-pi; >>x=fzero(fu,2) x= 2.4482 或 >>x=fsolve('0.8*x+atan(x)-pi',2) x= 2.4482 当然了,对于该方程也可以用第二种方法求解: >>x=solve('0.8*x+atan(x)-pi','x') x= 2.4482183943587910343011460497668 对于第一个例子,也可以用第三种方法求解: >>F=@(x)x^3-6*x^2-72*x-27 F= @(x)x^3-6*x^2-72*x-27 >>x=fzero(F,10) x= 12.1229 对于第二个例子,也可以用第三种方法: >>FUN=@(x)(x+2)^x-2 FUN= @(x)(x+2)^x-2 >>x=fzero(FUN,1) x= 0.6983 综上所述,可将常用的matlab解方程组的方法总结如下: 1、对于代数方程组Ax=b(A为系数矩阵,非奇异)的求解,MATLAB中有两种方法: (1)x=inv(A)*b — 采用求逆运算解方程组; (2)x=A\b — 采用左除运算解方程组。 例: x1+2x2=8 2x1+3x2=13 >>A=[1,2;2,3];b=[8;13]; >>x=inv(A)*b x = 2.00 3.00 >>x=A\b x = 2.00 3.00; 即二元一次方程组的解x1和x2分别是2和3。 2、用matlab解多次的方程组,有符号解法,方法是: 先解出符号解,然后用vpa(F,n)求出n位有效数字的数值解.具体步骤如下: 第一步: 定义变量syms x y z ...; 第二步: 求解[x,y,z,...]=solve('eqn1','eqn2',...,'eqnN','var1','var2',...'varN'); 第三步: 求出n位有效数字的数值解x=vpa(x,n);y=vpa(y,n);z=vpa(z,n);...。 如: 解二(多)元二(高)次方程组: x^2+3*y+1=0 y^2+4*x+1=0 解法如下: >>syms x y; >>[x,y]=solve('x^2+3*y+1=0','y^2+4*x+1=0'); >>x=vpa(x,4); >>y=vpa(y,4); 结果是: x = 1.635+3.029*i 1.635-3.029*i -.283 -2.987 y = 1.834-3.301*i 1.834+3.301*i -.3600 -3.307。 二元二次方程组,共4个实数根。 3、用matlab解高次方程组(非符号方程组)基本方法是: solve(s1,s2,…,sn,v1,v2,…,vn),即求表达式s1,s2,…,sn组成的方程组,求解变量分别v1,v2,…,vn。 具体例子如下: x^2 + x*y + y = 3 x^2 - 4*x + 3 = 0 解法: >> [x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3 = 0') 运行结果为 x = 1 3 y = 1 -3/2 即x等于1和3;y等于1和-1.5 或 >>[x,y] = solve('x^2 + x*y + y = 3','x^2 - 4*x + 3= 0','x','y') x = 1 3 y = 1 -3/2 结果一样,二元二方程都是4个实根。 通过这三个例子可以看出,用matlab解各类方程组都是可以的,方法也有多种,只是用到解方程组的函数,注意正确书写参数就可以了,非常方便。 第二部分解决实际问题的例子 【例一】MATLAB在结构化学中的应用举例 休克尔分子轨道的辅助计算: 休克尔分子轨道简称HMO,1931年由休克尔提出,是处理分子轨道以解决共轭分子的结构,探讨分子的性质和反应性能的半经验方法。 下面以丁二烯的HMO处理过程为例来说明MATLAB"求解矩阵的功能。 经休克尔基本假设化简后的久期方程为: =0 同除以b,令x= 得久期行列式 =0 应用MATLAB求解编程如下: >>symsx >>D=[x100;1x10;01x1;001x]; >>d=det(D) d=x^4-3*x^2+1 >>slove(‘x^4-3*x^2+1=0’) 得四解: 1/2*5^(1/2)+1/2 1/2-1/2*5^(1/2) 1/2*5^(1/2)-1/2 -1/2-1/2*5^(1/2) 由此便可得出分子轨道的能量等信息。 【例二】应用型实验: 工资问题 现有一个木工、一个电工、一个油漆工、一个装修工,共同合作完成他们各自的住房装修。 为了用工分配合理,决定每人工作13天,日工资不超过100元,且每人的总收入与总支出相等,试建立数学模型,以合理确定每人的日工资额。 工作天数分配方案 工种 天数 木工 电工 油漆工 装修工 在木工家的工作天数 2 1 6 3 在电工家的工作天数 4 5 2 4 在油漆工家的工作天数 4 4 3 2 在装修工家的工作天数 3 3 2 4 解: (1)问题分析与模型建立 木工、电工、油漆工、装修工,装修他们的房子,每人共工作13天,每人的日工资数应使得每人的总收入与总支出相等,这是一个投入产出问题。 设木工、电工、油漆工、装修工的日工资分别为: x1,x2,x3,x4。 则收支平衡关系式分别为: 2x1+x2+6x3+3x4=13x14x1+5x2+2x3+4x4=13x2 4x1+4x2+3x3+2x4=13x33x1+3x2+2x3+4x4=13x4 得数学模型: (2)Matlab求解 >>A=[-11163;2-412;22-51; 332-9]; >>B=rref(A) B= 1.000000-0.9912 01.00000-1.2719 001.0000-1.1053 0000 >>C=-B(: 4); >>C(4)=1; >>k=70; >>x=zeros(4,1); >>X=[]; >>whilemax(x)<=100 k=k+1; fori=1: 4 x(i)=C(i)*k; end X=[Xx]; end >>k=k-1 >>d=X(: k-70) 运行结果 k= 78 d= 77.3158 99.2105 86.2105 78.0000 齐次线性方程组的通解为: 由xi≤100得K=78,以确定木工、电工、油漆工 及装修工每人每天的日工资: x1=77x2=99x3=86x4=78 【例三】某地有三个产业,一个煤矿,一个发电厂和一条铁路,开采一元钱的煤,煤矿要支付0.25元的电费及0.25元的运输费;生产一元钱的电力,发电厂要支付0.65元的煤费,0.05元的电费及0.05元的运输费;创收一元钱的运输费,铁路要支付0.55元的煤费和0.10元的电费,在某一周内煤矿接到外地金额50000元定货,发电厂接到外地金额25000元定货,外界对地方铁路没有需求。 问三个企业间一周内总产值多少才能满足自身及外界需求? 三个企业间相互支付多少金额? 三个企业各创造多少新价值? 解: 这是一个投入产出分析问题。 设x1为本周内煤矿总产值,x2为电厂总产值,x3为铁路总产值,则 产出向量X= 外界需求向量D= 直接消耗矩阵C= 则原方程为(E-C)X=D 投入产出矩阵为B=C*diag(X) 总投入向量Y=ones(1,3)*B 新创造价值向量F=X-Y’ 消耗部门 外界需求 总产出 煤矿 电厂 铁路 生产部门 煤矿 0 36506 15582 50000 102088 电厂 25522 2808 2833 25000 56163 铁路 25522 2808 0 0 28330 新创造价值 51044 14041 9915 总产出 102088 56163 28330 【例四】在OH-离子作用下,硝基苯甲酸乙酯的纾解反应为: NO2C6H4COOC2H5+OH-=NO2C6H4COO-+C2H5OH 在15℃时,测得动力学数据如表1所示 表1硝基苯甲酸乙酯水解反应的转化率时间的变化(15℃) t/s α/% t/s α/% 120 32.95 330 58.05 180 41.75 530 69.00 240 48.80 600 70.35 α为酯水解转化率 解: 由二级反应速率微分公式 , 分离变量积分可得: , ~t为线性关系,可线性拟合,由斜率获得该二级反应的速率系数。 数据处理结果见表2。 表2硝基苯甲酸乙酯水解反应的 随时间的变化(15℃) t/s /dm3•mol-1 t/s /dm3•mol-1 120 9.828 330 27.68 180 14.34 530 44.51 240 19.06 600 47.45 运用Matlab编程线性拟合如下: clearall;clc t=[120180240330530600]; y=[9.82814.3419.0627.6844.5147.45]; p=polyfit(t,y,1) ti=linspace(t (1),t(end),100); yi=polyval(p,ti); plot(t,y,'o',ti,yi,'-') 计算结果如下: 斜率(k)=0.081267 保留4位有效数字,得反应速率系数k=0.08127dm3·mol-1·s-1。 【例五】: 在乙醇溶液中进行反应C2H5I+OH-=C2H5OH+I-,测得数据如表1所示,求反应的活化能。 表1不同温度下测得的反应速率系数 t/℃ k/10-3dm3·mol-1·s-1 t/℃ k/10-3dm3·mol-1·s-1 15.83 0.0503 59.75 6.71 32.02 0.368 90.61 119 解: 由阿仑尼乌斯微分方程dln{k}/dT=Ea/RT2,分离变量积分得: lnk=-EaR·T-1+lnk0。 Lnk~T-1为线性关系,可线性拟合,由斜率获得反应的活化能。 数据处理结果见表2: 表2按照阿伦尼乌斯公式处理得到的数据 ×103/K-1 ln(k/10-3dm3·mol-1·s-1) ×103/K-1 ln(k/10-3dm3·mol-1·s-1) 3.4606 -2.9898 3.0039 1.9036 3.2769 -0.9997 2.7491 4.7791 运用Matlab编程线性拟合如下: clearall;clc t=[3.46043.27693.00392.7491]; k=[0.05030.3686.71119]; k=log(k); p=polyfit(t,k,1) ti=linspace(t (1),t(end),100); ki=polyval(p,ti); plot(t,k,'o',ti,ki,'-') p*1000*8.3145 计算结果如下: 斜率(k)=-10892 保留4位有效数字,得Ea=-k·R=9.056×105J·mol-1。
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- Matlab 求解 线性方程组 非线性 方程组