Doolittle分解法matlab编程

Doolittle分解方法是一种用于解决线性方程组的数值方法，它可以将一个矩阵分解为一个下三角矩阵和一个上三角矩阵的乘积，从而可以方便地求解线性方程组。

在本文中，我们将介绍Doolittle分解方法的原理和实现过程，并用编程语言实现该方法来解方程组。

一、Doolittle分解方法原理1.1 Doolittle分解方法是一种LU分解的特例，它将一个矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积，即A=LU。

其中，L 的主对角线元素全为1，U的主对角线以上的元素全为0。

这样的分解可以方便地求解线性方程组Ax=b，其中b是一个已知的列向量。

1.2 Doolittle分解方法的具体实现过程是通过高斯消元法来实现的。

将矩阵A分解为一个下三角矩阵L和一个上三角矩阵U，然后通过回代法求解线性方程组Ax=b。

具体来说，我们首先将矩阵A分解为L 和U，然后用L和U的乘积代替原来的矩阵A，将原来的线性方程组Ax=b变为LUx=b，然后通过两次回代法求解线性方程组Ly=b和Ux=y，最终得到线性方程组的解x。

1.3 Doolittle分解方法的优点是可以方便地求解多个方程组，因为一旦矩阵A被分解为L和U，就可以通过多次回代法来求解不同的线性方程组，而不需要重新分解矩阵A。

1.4 Doolittle分解方法的缺点是需要对原始的矩阵A进行分解，这需要一定的计算量，特别是对于比较大的矩阵来说。

Doolittle分解方法在实际应用中往往需要结合其他数值方法来提高求解线性方程组的效率。

二、Doolittle分解方法的实现过程2.1 我们需要定义一个函数来实现Doolittle分解。

该函数的输入是一个矩阵A，输出是矩阵A的下三角矩阵L和上三角矩阵U。

2.2 接下来，我们需要通过高斯消元法来实现Doolittle分解。

具体来说，我们首先对矩阵A进行行变换和列变换，使得矩阵A的主对角线元素非零，然后逐步消去矩阵A的非主对角线元素，得到下三角矩阵L和上三角矩阵U。

第1章：绪论【1.2】设准确值为* 3.78694x =，*10y =，取它们的近似值分【1.1】按有效数字的定义，从两个方面说出1.0,1.00,1.000的不同含义【解】1.0,1.00,1.000的有效数字分别是两位，三位和四位；绝对误差限分别是0.05，0.005和0.0005别为123.7869, 3.780x x ==及129.9999, 10.1y y ==，试分析1212,,,x x y y 分别具有几位有效数字。

【解】*10.000040.00005x x -=<，1x 有5位有效数字；*20.006940.005x x -=>，2x 有2位有效数字；*10.000010.0005y y -=<，1y 有4位有效数字*2||0.10.5y y -=<，2y 有2位有效数字【1.3】（1）设p 的近似数有4位有效数字，求其相对误差限。

（2）用22/7和355/113作为 3.14159265p =L 的近似值，问它们各有几位有效数字？【解】（1）其绝对误差限是0.0005，则相对误差限为0.0005/3.1420.01591%r E ==（2）22/7 3.142857...=，有3位有效数字；355/113 3.14159292...=，有7位有效数字。

【1.4】试给出一种算法计算多项式32216180x a x a x a ++的函数值，使得运算次数尽可能少。

【解】24816328163281632012012,,,,x x x x x a x a x a x a x a x a x Þ++=++，总共8次乘法，两次加法【1.5】测量一木条长为542cm ，若其绝对误差不超过0.5cm ，问测量的相对误差是多少？【解】相对误差为0.5/5420.09%Î==【1.6】已知 2.71828e =L ，试问其近似值1232.7, 2.71, 2.718x x x ===各有几位有效数字？并给出他们的相对误差限。

《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名：学号：班级：2013年11月19日二阶椭圆偏微分方程第一边值问题摘要对于解二阶椭圆偏微分方程第一边值问题.课本上已经给出了相应的差分方程。

而留给我的难题就是把差分方程组表示成系数矩阵的形式.以及对系数进行赋值。

解决完这个问题之后.我在利用matlab 解线性方程组时.又出现“out of memory ”的问题。

因为99*99阶的矩阵太大.超出了分配给matlab 的使用内存。

退而求其次.当n=10.h=1/10或n=70.h=1/70时.我都得出了很好的计算结果。

然而在解线性方程组时.无论是LU 分解法或高斯消去法.还是gauseidel 迭代法.都能达到很高的精度。

关键字：二阶椭圆偏微分方程 差分方程 out of memory LU 分解 高斯消去法 gauseidel 迭代法一、题目重述解微分方程：()()2222((,))((,))()(,)()(,)(,)1y x x x y y x y yxxyxye u x y e u x y x y u x y x y u x y u x y y e x e e y x e--+++-+=-++++已知边界:(0,)1,(1,),(,0)1,(,1)y x u y u y e u x u x e ====求数值解, 把区域[0,1][0,1]G =?分成121/100,1/100h h ==.n =100 注：老师你给的题F 好像写错了.应该把22x y y e x e +改成22y x y e x e +。

二、问题分析与模型建立2.1微分方程上的符号说明()()22221y x xy xy y e x e e y x e -++++2.2课本上差分方程的缺陷课本上的差分方程为：举一个例子：当i=2,j=3时.；当i=3,j=3时.。

但是.显然这两个不是同一个数.其大小也不相等。

一、引言数值计算方法又称数值分析，是研究适合计算机求解的各种数学问题的近似方法及其理论。

它的内容包括函数逼近、数值微分与积分、非线性方程（组）的数值解、数值代数、常微分与偏微分方程数值解等。

这门课程起着承上启下的作用，承上是使线性代数、高等数学中的原理得以应用，方法得以实现，启下是为后续课程中数学问题的建模和求解提供方法，是高等理工科院校的重要基础课程。

如今，数值计算、理论研究及物理实验并列成为当今世界科学活动的三种主要方式。

为众多的科学与工程问题提供计算方法，提高计算的可靠性、有效性和精确性，是《数值计算方法》这门课程研究的主要内容。

在长期的教学实践中体会到在《数值计算方法》课程中做好理论内容的传授和学生实践能力的培养这两个环节非常重要。

如何合理的利用计算机软件进行有效地教与学是值得探讨的一个课题。

本文以具体教学为例，介绍了MATLAB软件在提高《数值计算方法》课堂教学质量中的具体使用。

二、MATLAB软件引入的必要性MATLAB是美国MathWorks公司自上世纪80年代中期推出的数学软件，其优秀的数值计算能力和卓越的数据可视化能力使其很快在数学软件中脱颖而出。

在欧美等高校，MATLAB已经成为线性代数、自动控制理论、数理统计、数值计算方法、动态系统仿真等高级课程的基本教学工具。

以前的《数值计算方法》课程常采用FORTRAN或者C语言进行教学和实验，要求学生既要对算法有充分的了解，又要熟练掌握这两种语言的语法和编程技巧，导致学生和教师将大量的时间和精力都花在烦琐的程序编写以及对各种结果的绘图上，学习效果往往令人不满意。

正如FORTRAN和C等高级语使人们摆脱了需要直接对计算机硬件资源进行操作一样，被称为第四代计算机语言的MATLAB，以其简洁的、更符合人们思维习惯的代码以及强大的绘图能力备受青睐。

《数值计算方法》课程内容多、课时少，如果运用传统教学方法，有些内容得不到细致地讲解，易使学生产生厌学情绪，收不到良好的教学效果。

matlab doolittle分解法代码在MATLAB中，Doolittle分解法（也称为LU分解）是一种常用的求解线性方程组的方法。

以下是一个使用MATLAB实现Doolittle分解法的示例代码：```MATLAB生成一个随机矩阵AA = rand(3, 3);初始化L和U矩阵为0L = eye(3);U = eye(3);采用Doolittle分解法分解矩阵Afor i = 1:3for j = 1:3if i == jU(i, j) = A(i, j);elseL(i, j) = A(i, j);endend寻找最大值并交换行max_idx = find(max(A(:, i)));L(:, i) = L(:, i) / A(max_idx, i);U(:, i) = U(:, i) / A(max_idx, i);end打印L和U矩阵disp('L矩阵：');disp(L);disp('U矩阵：');disp(U);验证分解是否正确A_original = A;L_original = eye(3);U_original = eye(3);[L_original, U_original] = doolittle(A_original);打印原始的L和U矩阵disp('原始的L矩阵：');disp(L_original);disp('原始的U矩阵：');disp(U_original);计算分解后的矩阵乘积与原矩阵是否相等result = L_original * U_original;disp('分解后的矩阵：');disp(result);检查分解后的矩阵是否与原矩阵相等isequal(A_original, result)```以上代码首先生成一个3x3的随机矩阵A，然后使用Doolittle分解法对其进行分解，并打印出分解后的L和U矩阵。

数值计算⽅法matlab程序function [x0,k]=bisect1(fun1,a,b,ep)if nargin<4ep=1e-5;endfa=feval(fun1,a);fb=feval(fun1,b);if fa*fb>0x0=[fa,fb];k=0;return;endk=1;while abs(b-a)/2>epx=(a+b)/2;fx=feval(fun1,x);if fx*fa<0b=x;fb=fx;elsea=x;fa=fx;k=k+1;endendx0=(a+b)/2;>> fun1=inline('x^3-x-1');>> [x0,k]=bisect1(fun1,1.3,1.4,1e-4)x0 =1.3247k =7>>N=500;endif nargin<3ep=1e-5;endx=x0;x0=x+2*ep;while abs(x-x0)>ep & kx0=x;x=feval(fun1,x0);k=k+1;endx0=x;if k==Nwarning('已达最⼤迭代次数')end>> fun1=inline('(x+1)^(1/3)');>> [x0,k]=iterate1(fun1,1.5)x0 =1.3247k =7>> fun1=inline('x^3-1');>> [x0,k]=iterate1(fun1,1.5)x0 =Infk =9>>Steffesen加速迭代（简单迭代法的加速）function [x0,k]=steffesen1(fun1,x0,ep,N) if nargin<4N=500;endx=x0;x0=x+2*ep;k=0;while abs(x-x0)>ep & kx0=x;y=feval(fun1,x0);z=feval(fun1,y);x=x0-(y-x0)^2/(z-2*y+x0);k=k+1;endx0=x;if k==Nwarning('已达最⼤迭代次数')end>> fun1=inline('(x+1)^(1/3)');>> [x0,k]=steffesen1(fun1,1.5)x0 =1.3247k =3>> fun1=inline('x^3-1');>> [x0,k]=steffesen1(fun1,1.5)x0 =1.3247k =6Newton迭代function [x0,k]=Newton7(fname,dfname,x0,ep,N) if nargin<5N=500;endendx=x0;x0=x+2*ep;k=0;while abs(x-x0)>ep & kx0=x;x=x0-feval(fname,x0)/feval(dfname,x0);k=k+1;endx0=x;if k==Nwarning('已达最⼤迭代次数')end>> fname=inline('x-cos(x)');>> dfname=inline('1+sin(x)');>> [x0,k]=Newton7(fname,dfname,pi/4,1e-8) x0 =0.7391k =4⾮线性⽅程求根的Matlab函数调⽤举例：1.求多项式的根：求f(x)=x^3-x-1=0的根：>> roots([1 0 -1 -1])ans =1.3247-0.6624 + 0.5623i-0.6624 - 0.5623i2.求⼀般函数的根>> fun=inline('x*sin(x^2-x-1)','x')fun =Inline function:fun(x) = x*sin(x^2-x-1)>> fplot(fun,[-2 0.1]);grid on-1.5956>> x=fzero(fun,[-1 -0.1])x =-0.6180[x,f,h]=fsolve(fun,-1.6)x =-1.5956f =1.4909e-009h =1（h>0表⽰收敛，h<0表⽰发散，h=0表⽰已达到设定的计算函数值的最⼤次数）第三章：线性⽅程组的数值解法1. ⾼斯消元法function [A,x]=gauss3(A,b)%本算法⽤顺序⾼斯消元法求解线性⽅程组n=length(b);A=[A,b];for k=1:n-1A((k+1):n,(k+1):(n+1))=A((k+1):n,(k+1):(n+1))-A((k+1):n,k)/A(k,k)*A(k,(k+1):(n+1)); A((k+1):n,k)=zeros(n-k,1);A;endx=zeros(n,1);%上⾯为消元过程x(n)=A(n,n+1)/A(n,n);for k=n-1:-1:1x(k)=(A(k,n+1)-A(k,(k+1):n)*x((k+1:n)))/A(k,k);end%上⾯为回代过程>> A=[2 3 4;3 5 2;4 3 30];>> b=[6,5,32]'b =>> [A,x]=gauss3(A,b)A =2.00003.00004.0000 6.00000 0.5000 -4.0000 -4.00000 0 -2.0000 -4.0000x =-1382列选主元的⾼斯消元法：function [A,x]=gauss5(A,b)%本算法⽤列选主元的⾼斯消元法求解线性⽅程组n=length(b);A=[A,b];for k=1:n-1%选主元[ap,p]=max(abs(A(k:n,k)));p=p+k-1;if p>kt=A(k,:);A(k,:)=A(p,:);A(p,:)=t;end%消元A((k+1):n,(k+1):(n+1))=A((k+1):n,(k+1):(n+1))-A((k+1):n,k)/A(k,k)*A(k,(k+1):(n+1)); A((k+1):n,k)=zeros(n-k,1);end%回代x=zeros(n,1);x(n)=A(n,n+1)/A(n,n);>> A=[2 3 4;3 5 2;4 3 30]; b=[6,5,32]';>> [A,x]=gauss5(A,b)A =4.0000 3.0000 30.0000 32.00000 2.7500 -20.5000 -19.00000 0 0.1818 0.3636x =-1382三⾓分解法：Doolittle 分解function [L,U]=doolittle1(A)n=length(A);U=zeros(n);L=eye(n);U(1,:)=A(1,:);L(2:n,1)=A(2:n,1)/U(1,1);for k=2:nU(k,k:n)=A(k,k:n)-L(k,1:k-1)*U(1:k-1,k:n);L(k+1:n,k)=A(k+1:n,k)-L(k+1:n,1:k-1)*U(1:k-1,n)/U(k,k); End y=zeros(n,1);x=y;y(1)=b(1);for i=2:ny(i)=b(i)-L(i,1:i-1)*y(1:i-1);endx(n)=y(n)/U(n,n);for i=n-1:-1:1x(i)=(y(i)-U(i,i+1:n)*x(i+1:n))/U(i,i);end>> A=[1 2 3;2 5 2 ;3 1 5];b=[14 18 20]';>> [L,U,x]=doolittle1(A,b)3 -8 1U =1 2 30 1 -40 0 -36x =2.83331.33332.8333平⽅根法：function [L,x]=choesky3(A,b)n=length(A);L=zeros(n);L(:,1)=A(:,1)/sqrt(A(1,1));for k=2:nL(k,k)=A(k,k)-L(k,1:k-1)*L(k,1:k-1)';L(k,k)=sqrt(L(k,k));for i=k+1:nL(i,k)=(A(i,k)-L(i,1:k-1)*L(k,1:k-1)')/L(k,k); endendy=zeros(n,1);x=y;y(1)=b(1)/L(1,1);for i=2:ny(i)=(b(i)-L(i,1:i-1)*y(1:i-1))/L(i,i);endx(n)=y(n)/L(n,n);for i=n-1:-1:1x(i)=(y(i)-L(i+1:n,i)'*x(i+1:n))/L(i,i);end>> A=[4 -1 1;-1 4.25 2.75;1 2.75 3.5]-1.0000 4.2500 2.75001.00002.75003.5000>> b=[4 6 7.25]'b =4.00006.00007.2500[L,x]=choesky3(A,b)L =2.0000 0 0-0.5000 2.0000 00.5000 1.5000 1.0000x =111>>迭代法求⽅程组的解Jacobi迭代法：function [x,k]=jacobi2(a,b,x0,ep,N)%本算法⽤Jacobi迭代求解ax=b,⽤分量形式n=length(b); k=0;if nargin<5N=500;endif nargin<4ep=1e-5;endif nargin<3x0=zeros(n,1);y=zeros(n,1);while norm(x-x0,inf)>ep & kk=k+1;x0=x;for i=1:ny(i)=b(i);for j=1:nif j~=iy(i)=y(i)-a(i,j)*x0(j);endendif abs(a(i,i))<1e-10|k==Nwarning('a(i,i) is too small');returnendy(i)=y(i)/a(i,i);endx=y;enda=[4 3 0;3 4 -1; 0 -1 4];b=[24 30 -24]'；[x,k]=jacobi2(a,b)x =3.00004.0000-5.0000k =59Gauss-seidel迭代法：function [x,k]=gaussseide2(a,b,x0,ep,N)%本算法⽤Gauss-seidel迭代求解ax=b,⽤分量形式n=length(b); k=0;if nargin<5N=500;endendif nargin<3x0=zeros(n,1);y=zeros(n,1);endx=x0;x0=x+2*ep;while norm(x-x0,inf)>ep & kk=k+1;x0=x;y=x;for i=1:nz(i)=b(i);for j=1:nif j~=iz(i)=z(i)-a(i,j)*x(j);endendif abs(a(i,i))<1e-10|k==Nwarning('a(i,i) is too small');returnendz(i)=z(i)/a(i,i);x(i)=z(i);endend[x,k]=gaussseide2(a,b)x =3.00004.0000-5.0000k =25最速下降法function [x,k]=zuisuxiajiang(A,b,x0,ep,N)N=500;endif nargin<4ep=1e-8;endif nargin<3x0=ones(n,1);endx=x0;x0=x+2*ep;r=b-A*x;d=r;k=0;while norm(x-x0,inf)>ep & kk=k+1;x0=x;lamda=(d'*d)/(d'*A*d);x=x0+lamda*d;r=b-A*x;d=r;endif k==Nwarning('已达最⼤迭代次数')end共轭梯度算法function [x,k]=gongertidufa(A,b,x0,ep,N) %本算法⽤共轭梯度算法求解正定⽅程组Ax=b,,n=length(b);if nargin<5N=500;endif nargin<4ep=1e-8;x0=x+2*ep;r=b-A*x;d=r;k=0;while norm(x-x0,inf)>ep & kx0=x;lamda=(r'*r)/(d'*A*d);r1=r;x=x0+lamda*d;r=b-A*x;beta=(r'*r)/(r1'*r1);d=r+beta*d;endif k==Nwarning('已达最⼤迭代次数') end常微分⽅程数值解function [x,y]=Euler1(fun,xspan,y0,h)%本算法⽤欧拉格式计算微分⽅程y'=f(x,y)的解。

1．列主元高斯消去法M文件function[x]=gauss(a,b)n=length(a);x=zeros(n,1);a=[a b];for k=1:n-1max=k;for i=k+1:nif a(i,k)>a(max,k)max=i;endendtemp=a(k,k:n+1);a(k,k:n+1)=a(max,k:n+1);a(max,k:n+1)=temp;for i=k+1:na(i,k)=-a(i,k)/a(k,k);a(i,k+1:n+1)=a(i,k+1:n+1)+a(i,k)*a(k,k+1:n+1);endendx(n,1)=a(n,n+1)/a(n,n);for i=n-1:-1:1sum=0;for j=i+1:nsum=sum+x(j,1)*a(i,j);endx(i,1)=(a(i,n+1)-sum)/a(i,i);endMatlab运行结果2．LU三角分解法M文件function y=LU(A,B);n=length(A);A=[A B];for k=1:n-1;for i=k:n;if(abs(A(i,k))==max(abs(A(k:n,k))))P(k)=i;temp=A(k,:);A(k,:)=A(i,:);A(i,:)=temp;endendfor j=k+1:n;A(j,k)=A(j,k)/A(k,k);A(j,k+1:n+1)=A(j,k+1:n+1)-A(j,k)*A(k,k+1:n+1);endendP(n)=n;L(1,1)=1;L(2:n,1)=A(2:n,1);L(1,2:n)=0;U(1,1)=A(1,1);U(2:n,1)=0;U(1,2:n)=A(1,2:n);for i=2:n;L(i,1:i-1)=A(i,1:i-1);L(i,i)=1;L(i,i+1:n)=0;U(i,1:i-1)=0;U(i,i:n)=A(i,i:n);endx(n) = A(n,n+1)/U(n,n);for k = n-1:-1:1x(k)=A(k,n+1);for p=n:-1:k+1;x(k) = x(k)-U(k,p)*x(p); endx(k)=x(k)/U(k,k);endxLUPEndMatlab运行结果3．龙贝格(Romberg)算法M文件function[t]=romberg(f,a,b,e)t=zeros(15,4);t(1,1)=(b-a)/2*(f(a)+f(b));for k=2:4sum=0;for i=1:2^(k-2)sum=sum+f(a+(2*i-1)*(b-a)/2^(k-1));endt(k,1)=0.5*t(k-1,1)+(b-a)/2^(k-1)*sum;for i=2:kt(k,i)=(4^(i-1)*t(k,i-1)-t(k-1,i-1))/(4^(i-1)-1);endendfor k=5:15sum=0;for i=1:2^(k-2)sum=sum+f(a+(2*i-1)*(b-a)/2^(k-1));endt(k,1)=0.5*t(k-1,1)+(b-a)/2^(k-1)*sum;for i=2:4t(k,i)=(4^(i-1)*t(k,i-1)-t(k-1,i-1))/(4^(i-1)-1);endif k>6if abs(t(k,4)-t(k-1,4))<edisp(['答案',num2str(t(k,4))]);break;endendendif k>=15disp(['溢出']);endMatlab运行结果4．最小二乘法M文件function[a,max,det]=zuixiaoerchengfa(x,y,r) n=length(x);c=ones(n,r+1);for i=2:r+1for j=1:nc(j,i)=x(j)^(i-1);endendA=c'*c;b=c'*y';a=inv(A)*b;det=0;max=0;for i=1:nsum=a(1);for j=2:r+1sum=sum+a(j)*x(i)^(j-1);endcc=abs(y(i)-sum);if cc>maxmax=cc;enddet=det+cc^2;enddet=sqrt(det);Matlab运行结果§2.1.1 二分法的MATLAB主程序function [k,x,wuca,yx]=erfen(a,b,abtol)a(1)=a; b(1)=b;ya=fun(a(1)); yb=fun(b(1)); %程序中调用的fun.m 为函数if ya* yb>0,disp('注意：ya*yb>0,请重新调整区间端点a和b.'), returnendmax1=-1+ceil((log(b-a)- log(abtol))/ log(2)); % ceil是向+方向取整∞for k=1: max1+1a;ya=fun(a); b;yb=fun(b); x=(a+b)/2;yx=fun(x); wuca=abs(b-a)/2; k=k-1;[k,a,b,x,wuca,ya,yb,yx]if yx==0a=x; b=x;elseif yb*yx>0b=x;yb=yx;elsea=x; ya=yx;endif b-a< abtol , return, endendk=max1; x; wuca; yx=fun(x);§2.1.2 不动点迭代法的MATLAB主程序function [k,piancha,xdpiancha,xk,yk]=diedai2(x0,tol,ddmax) x(1)=x0;for i=1: ddmaxx(i+1)=fun(x(i));piancha=abs(x(i+1)-x(i));xdpiancha=piancha/( abs(x(i+1))+eps);i=i+1;xk=x(i);yk=fun(x(i)); [(i-1) piancha xdpiancha xk yk]if (piancha<tol)|(xdpiancha< tol)k=i-1; xk=x(i);return;endendif i>ddmaxdisp('迭代次数超过给定的最大值ddmax')k=i-1; xk=x(i);yk=fun(x(i));[(i-1) piancha xdpiancha xk yk];return;endP=[(i-1),piancha,xdpiancha,xk,yk]';§2.1.3 艾特肯加速迭代法的MATLAB主程序function [k,xk,yk,p]= Aitken (x0,tol, ddmax)x(1)=x0;for i=1: ddmaxx1(i+1)=fun(x(i)); x2(i+1)=fun(x1(i+1));x(i+1)=x2(i+1)-(x2(i+1)-x1(i+1))^2/(x2(i+1)-2*x1(i+1)+ x(i));piancha=abs(x(i+1)-x(i));xdpiancha= piancha/( abs(x(i+1))+eps);i=i+1; xk=x(i);yk=fun(x(i));if (piancha<tol)|(xdpiancha<tol)k=i-1; xk=x(i); yk=fun(x(i));m=[0,1:i-1]; p=[m',x1',x2',x'];return;endendif i>ddmaxdisp('迭代次数超过给定的最大值ddmax')k=i-1; xk=x(i); yk=fun(x(i));m=[0,1:i-1]; p=[m',x1',x2',x'];return;endm=[0,1:i-1]; p=[m',x1',x2',x'];§2.1.4 牛顿切线法的MATLAB主程序function [k,xk,yk,piancha,xdpiancha]=newtonqx(x0,tol,ftol,gxmax) x(1)=x0;for i=1: gxmaxx(i+1)=x(i)-fun(x(i))/(dfun(x(i))+eps);piancha=abs(x(i+1)-x(i));xdpiancha= piancha/( abs(x(i+1))+eps); i=i+1;xk=x(i);yk=fun(x(i)); [(i-1) xk yk piancha xdpiancha]if (abs(yk)<ftol)&((piancha<tol)|(xdpiancha< tol))k=i-1; xk=x(i);[(i-1) xk yk piancha xdpiancha]return;endendif i>gxmaxdisp('请注意：迭代次数超过给定的最大值gxmax。

科学计算—理论、方法及其基于ＭＡＴＬＡＢ的程序实现与分析 三、 解线性方程组（线性矩阵方程）解线性方程组是科学计算中最常见的问题。

所说的“最常见”有两方面的含义：１） 问题的本身是求解线性方程组；２） 许多问题的求解需要或归结为线性方程组的求解。

关于线性方程组B A x B Ax 1-=⇒=（１）其求解方法有两类：１） 直接法：高斯消去法（Gaussian Elimination ）； ２） 间接法：各种迭代法（Iteration ）。

１、高斯消去法１） 引例考虑如下（梯形）线性方程组：()⎪⎩⎪⎨⎧==+==+-=⇒⎪⎪⎪⎭⎫ ⎝⎛=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--⇔⎪⎩⎪⎨⎧==-=+-5.0141315.3221122004301211214322332321321332321x x x x x x x x x x x x x x x 高斯消去法的求解思路：把一般的线性方程组（１）化成（上或下）梯形的形式。

２）高斯消去法——示例考虑如下线性方程组：⎪⎪⎪⎭⎫ ⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---⇔⎪⎩⎪⎨⎧=++-=-+-=+-306015129101.2001.221113060129501.2001.221321321321321x x x x x x x x x x x x １） 第一个方程的两端乘12加到第二个方程的两端，第一个方程的两端乘－１加到第三个方程的两端，得⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--3060031110001.0001.00111321x x x２） 第二个方程的两端乘001.010-加到第三个方程的两端，得 ⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡--60600311010001.0001.00111321x x x３） 从上述方程组的第三个方程依此求解，得()⎪⎩⎪⎨⎧==+-==+-=600300001.03100024011332321x x x x x x ３）高斯消去法的不足及其改进——高斯（全、列）主元素消去法在上例中，由于建模、计算等原因，系数2.001而产生0.0005的误差，实际求解的方程组为⎪⎪⎪⎭⎫⎝⎛-=⎪⎪⎪⎭⎫ ⎝⎛⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡---306015129101.20005.22111321x x x ⎪⎩⎪⎨⎧===⇒70.4509.30142.2565321x x x注：数值稳定的算法高斯列主元素消去法就是在消元的每一步选取（列）主元素—一列中绝对值最大的元取做主元素，高斯列主元素消去法是数值稳定的方法。