Matlab求解线性方程组、非线性方程组.docx
- 格式:docx
- 大小:10.70 KB
- 文档页数:6
本科毕业论文(2010 届)题目线性方程组的直接解法及matlab的实现学院数学与信息工程学院专业数学与应用数学班级2006级数学1 班学号0604010127学生姓名胡婷婷指导教师王洁完成日期2010年5月摘要随着科技技术的发展及人类对自然界的不断探索模拟。
在自然科学和工程问题中的很多问题的解决常常归结为线性代数问题!本文的主要内容是对线性方程组求解方法的探讨,主要介绍了四种求解线性方程组的方法,第一种是教科书上常见的消元法,我们称之为基本法。
第二种方法是标准上三角形求解法,即将增广矩阵经过初等变换后化成标准上三角形,然后求解.它改进了一般教科书上的常见方法,与常见方法比较有如下优点:1)规范了自由未知量的选取;2)只用矩阵运算;3)减少了计算量.第三种方法是对特定的方程组(系数矩阵A为n阶对称正定矩阵,且A的顺序主子式均不为零。
)的求解方法进行描述,并且为这种线性方程的求解提供了固定的公式化的方法。
第四种方法是对现在实际问题中常常会遇到的系数矩阵为三对角矩阵的方程组的求解方法。
同时给出这几种方法的数值解法(matlab程序),由于运用电脑软件求解,所以必须考虑计算方法的时间、空间上的效率以及算法的数值稳定性问题,所以针对不同类型的线性方程组有不同的解法.但是,基本的方法可以归结为两大类,即直接法和迭代法.关键词高斯消去法;三角分解法;乔莱斯基分解法;追赶法AbstractSystems of linear equations are associated with many problems in engineering and scinence ,as well as with applications of mathematics to the social sciences and the quantitative study of business and economic problems.The main content of this article is the method for solving linear equations,we introduce four methods for solving linear equations in this paper。
求解线性方程组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的解;X=B/A表示矩阵方程XA=B的解。
对方程组X=A\B,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A的列数,方程X=B/A同理。
如果矩阵A不是方阵,其维数是m×n,则有:m=n 恰定方程,求解精确解;m>n 超定方程,寻求最小二乘解;m<n 不定方程,寻求基本解,其中至多有m个非零元素。
实验2 利用matlab 解(非)线性、微分方程(组)
一、实验目的
1、线性方程组的解法:直接求解法和迭代法;
2、非线性方程以及非线性方程组的求解;
3、微分方程的数值解。
二、实验内容
1、对于下列线性方程组:
⎥⎥⎥⎦
⎤⎢⎢⎢⎣⎡=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡66136221143092x
(1) 请用直接法求解;
(2) 请用LU 分解方法求解;
(3) 请用QR 分解方法求解;
(4) 请用Cholesky 分解方法求解。
2、设迭代精度为10-6,分别用Jacobi 迭代法、Gauss-Serdel 迭代法求解下列线性方程组,并比较此两种迭代法的收敛速度。
⎪⎩
⎪⎨⎧=+-=-+-=-510272109103232121x x x x x x x 3、求解非线性方程010=-+-x xe x 在2附近的根。
4、求下列非线性方程组在(0.5,0.5) 附近的数值解。
⎪⎩
⎪⎨⎧=-+=-+02)sin(02)cos(y x xe y ye x 5、通过画图方法描述某非刚性体的运动方程的微分方程⎪⎩⎪⎨⎧=-==2
1131232151.0y y y y y y
y y y ,其初始条件为⎪⎩⎪⎨⎧===1)0(1)0(0)0(3
21y y y 。
6、求二阶微分方程)2sin(3t y e y t y t =-'+'', 1)0(=y ,1)0(-='y 在20≤≤t 时的数值图解。
m a t l a b非线性方程求解要点非线性方程的解法1引 言数学物理中的许多问题归结为解函数方程的问题,即,0)(=x f (1.1)这里,)(x f 可以是代数多项式,也可以是超越函数。
若有数*x 为方程0)(=x f 的根,或称函数)(x f 的零点。
设函数)(x f 在],[b a 内连续,且0)()(<b f a f 。
根据连续函数的性质知道,方程0)(=x f 在区间],[b a 内至少有一个实根;我们又知道,方程0)(=x f 的根,除了极少简单方程的根可以用解析式表达外,一般方程的根很难用一个式子表达。
即使能表示成解析式的,往往也很复杂,不便计算。
所以,具体求根时,一般先寻求根的某一个初始近似值,然后再将初始近似值逐步加工成满足精度要求为止。
如何寻求根的初始值呢?简单述之,为了明确起见,不妨设)(x f 在区间],[b a 内有一个实的单根,且0)(,0)(><b f a f 。
我们从左端出点a x =0出发,按某一预定的步长h 一步一步地向右跨,每跨一步进行一次根的“搜索”,即检查每一步的起点k x 和1+k x (即,h x k +)的函数值是否同号。
若有:0)(*)(≤+h x f x f k k (1.2)那么所求的根必在),(h x x k k +内,这时可取k x 或h x k +作为根的初始近似值。
这种方法通常称为“定步长搜索法”。
另外,还是图解法、近似方程法和解析法。
2 迭代法2.1 迭代法的一般概念迭代法是数值计算中一类典型方法,不仅用于方程求根,而且用于方程组求解,矩阵求特征值等方面。
迭代法的基本思想是一种逐次逼近的方法。
首先取一个精糙的近似值,然后用同一个递推公式,反复校正这个初值,直到满足预先给定的精度要求为止。
对于迭代法,一般需要讨论的基本问题是:迭代法的构造、迭代序列的收敛性天收敛速度以及误差估计。
这里,主要看看解方程迭代式的构造。
第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 方程组解一、概述Matlab是一种强大的数学计算软件,它可以用来解决各种数学问题,包括解方程组。
在Matlab中,求解方程组是一个非常重要的功能,因为很多实际问题都可以转化为方程组的形式。
本文将详细介绍如何使用Matlab求解线性方程组和非线性方程组。
二、线性方程组1. 线性方程组的定义线性方程组是指各个未知量的次数都不超过1次的代数方程组。
例如:2x + 3y = 54x - 5y = 6就是一个包含两个未知量x和y的线性方程组。
2. Matlab中求解线性方程组方法在Matlab中,可以使用“\”或者“inv()”函数来求解线性方程组。
其中,“\”表示矩阵左除,即Ax=b时,求解x=A\b;“inv()”函数表示矩阵求逆,即Ax=b时,求解x=inv(A)*b。
例如,在Matlab中求解以下线性方程组:2x + 3y = 54x - 5y = 6可以使用以下代码:A=[2,3;4,-5];b=[5;6];x=A\b输出结果为:x =1.00001.0000其中,“A”为系数矩阵,“b”为常数矩阵,“x”为未知量的解。
三、非线性方程组1. 非线性方程组的定义非线性方程组是指各个未知量的次数超过1次或者存在乘积项、幂项等非线性因素的代数方程组。
例如:x^2 + y^2 = 25x*y - 3 = 0就是一个包含两个未知量x和y的非线性方程组。
2. Matlab中求解非线性方程组方法在Matlab中,可以使用“fsolve()”函数来求解非线性方程组。
该函数需要输入一个函数句柄和初始值向量,输出未知量的解向量。
例如,在Matlab中求解以下非线性方程组:x^2 + y^2 = 25x*y - 3 = 0可以使用以下代码:fun=@(x)[x(1)^2+x(2)^2-25;x(1)*x(2)-3];x0=[1;1];[x,fval]=fsolve(fun,x0)输出结果为:Local minimum found.Optimization completed because the size of the gradient is less thanthe default value of the function tolerance.<stopping criteria details>ans =1.60561.8708其中,“fun”为函数句柄,表示要求解的非线性方程组,“x0”为初始值向量,“[x,fval]”为输出结果,其中“x”表示未知量的解向量,“fval”为函数值。
用matlab进行方程组求解的案例一、背景介绍方程组求解是数学中一个重要的问题,也是许多科学和工程领域中常见的问题。
在过去,人们通常使用手算或者计算器进行方程组求解,但是这种方法效率低下且容易出错。
随着计算机的发展和普及,使用计算机进行方程组求解已经成为了一种主流的方法。
Matlab是一款强大的数值计算软件,它提供了许多用于求解方程组的函数和工具箱。
本文将介绍如何使用Matlab进行方程组求解,并通过一个实例来说明其具体应用。
二、Matlab中的方程组求解函数在Matlab中,有多种函数可以用于求解方程组,包括线性方程组、非线性方程组、常微分方程等。
下面列出一些常用的函数:1. linsolve:用于求解线性方程组;2. fsolve:用于求解非线性方程组;3. ode45:用于求解常微分方程。
三、实例介绍假设有如下非线性方程组:x^2 + y^2 = 1x^3 - y = 0我们需要使用Matlab对其进行求解。
四、代码实现1. 定义函数首先需要定义一个函数,输入为变量向量x=[x,y],输出为方程组的值向量f=[f1,f2]。
代码如下:function f = myfun(x)f(1) = x(1)^2 + x(2)^2 - 1;f(2) = x(1)^3 - x(2);end2. 求解方程组使用fsolve函数求解方程组。
代码如下:x0 = [0,0]; % 初始值[x,fval] = fsolve(@myfun,x0);其中,x为方程组的解向量,fval为方程组的值向量。
3. 结果展示将求解结果输出。
代码如下:fprintf('x=%f, y=%f\n',x(1),x(2));运行程序后,得到如下结果:x=0.682327, y=0.731689五、结论本文介绍了如何使用Matlab进行非线性方程组求解,并通过一个实例进行了说明。
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;。
MATLAB教程方程组MATLAB(Matrix Laboratory)是一种高级计算机语言和环境,被广泛应用于科学、工程和其他技术领域。
它具有强大的矩阵操作能力,以及丰富的函数库,能够快速、高效地完成各种计算任务。
在MATLAB中,方程组的求解是一项重要的任务,本教程将介绍如何使用MATLAB求解线性和非线性方程组。
1.线性方程组的求解线性方程组是指方程中的未知量只有一次出现,并且未知量之间的关系是线性的。
在MATLAB中,可以使用“\”运算符求解线性方程组。
首先,定义一个线性方程组。
例如,我们要求解以下方程组:2x+4y+z=103x+2y-z=-7x-y+2z=5可以将方程组的系数矩阵和常数矩阵分别定义为A和B:A=[2,4,1;3,2,-1;1,-1,2];B=[10;-7;5];然后,使用“\”运算符求解方程组,并将结果赋值给未知量向量X:X=A\B;最后,打印出未知量向量X的值:disp(X);这样,就可以得到方程组的解。
在上述例子中,解为X=[1;-2;3]。
2.非线性方程组的求解非线性方程组是指方程中的未知量出现在非线性函数中,未知量之间存在复杂的关系。
在MATLAB中,可以使用fsolve函数求解非线性方程组。
首先,定义一个非线性方程组。
例如,我们要求解以下方程组:x^2+y^2=25x * cos(y) + y * sin(x) = 10然后,定义一个匿名函数,将方程组的函数表达式作为输入参数:接下来,使用fsolve函数求解方程组,并将初始猜测值[1, 1]作为输入参数:initialGuess = [1; 1];solution = fsolve(equations, initialGuess);最后,打印出方程组的解:disp(solution);这样,就可以得到非线性方程组的解。
在上述例子中,解为solution = [3.9781; 2.3382]。
除了fsolve函数外,MATLAB还提供了其他用于求解非线性方程组的函数,如lsqnonlin、fsolve、fmincon等。
求解线性方程组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的解;
X= B/A表示矩阵方程XA=B的解。
对方程组X = A∖B ,要求A和B用相同的行数,X和B有相同的列数,它的行数等于矩阵A 的列数,方程X= B/A 同理。
如果矩阵A不是方阵,其维数是m× n,则有: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 变换的基础上,其算法可靠性稍逊与奇异值求解,但速度较快;
【例7】
求解超定方程组
A=[2 -1 3;3 1 -5;4 -1 1;1 3 -13]
A=
2 -1 3
3 1 -5
4 -1 1
1 3 -13
b = [3 0 3 -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
可见x1 并不是方程Ax=b 的精确解,用x2=pinv(A)*b 所得的解与x1 相同。
三.欠定方程组欠定方程组未知量个数多于方程个数,但理论上有无穷个解。
MATLAB 将寻求一个基本解,其中最多只能有m 个非零元素。
特解由列主元qr 分解求得。
【例8】
解欠定方程组
A = [1-2 1 1;1 -2 1 -1;1 -2 1 5]
A=
1 -
2 1 1
1 -
2 1 -1
1 -
2 1 -1
1 -
2 1 5
b=[1 -1 5] '
x1=A\b
Warning:Rank deficient,rank=2 tol=4.6151e-015 x1= 0
-0.0000
1.0000
x2=pinv(A)*b
x2=
-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)>O 时,w(i)O,同时返回一个双向量w。
【例9】求方程组的非负最小二乘解
A=[3.4336 -0.5238 0.6710
-0.5238 3.2833 -0.7302
0.6710 -0.7302 4.0261];
b=[-1.000 1.5000 2.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。