下载文档原格式
一. 实验目的学习使用Guass 消去法和迭代法求解线性方程组二. 算法介绍或方法基础§3.1 Jacobi 迭代公式设方程组AX=b , 通过分离变量的过程建立Jacobi 迭代公式，即11,,0(1,2,,)1()(1,2,,)nijj i ii j ni i ij j j j iiiax b a i n x b a x i n a ==≠=≠==-=∑∑由此我们可以得到Jacobi 迭代公式：(1)1,1()(1,2,,)nk k i i ij j j j iii x b a x i n a +=≠=-=∑ §3.2 高斯-塞德尔迭代法(AX=b)(1)(1)(1)121,,,k k k i x x x +++- ()()()121,,,k k k i x x x - 的值，而Jacobi 迭代公式并不利用这些最新的近似值计算，仍用(1(=11)(1)11)(1,2,,)i nk k k ii ij j ij j j j i ii xb a x a x i n a -++==+--=∑∑ 这启发我们可以对其加以改进,即在每个分量的计算中尽量利用最新的迭代值，得到上式称为Gauss-Seidel 迭代法.1.矩阵的乘法>> a=[1 2;3 4];>> b=[5 5;7 8];>> a^2*b2.矩阵的除法>> a=[1 2 3;4 5 6;7 8 9];>> b=[1 0 0;0 2 0;0 0 3];>> a\b>> a/b3.使用冒号选出指定的元素>> a=[1 2 3;4 5 6;7 8 9];>> a(1:2,3)>> a(2:3,:)4. 线性方程组数值求解（1）Guass消去法clear;clc;%这个程序可适用于求解任何符合条件n元一次方程组A=input('请输入系数方阵A=');[m,n]=size(A);while m~=n | det(A)==0A=input('请输入系数方阵且det(A)~=0 A='); endB=input('请输入列向量B=');c=[A B];disp('合并后的增广矩阵为');cdisp('开始进行列主元消元法');for i=1:m-1for k=i+1:mif c(k,i)>c(i,i)|-c(k,i)>c(i,i)b=c(k,:);c(k,:)=c(i,:);c(i,:)=b;endcfor j=i+1:mc(j,:)=c(j,:)-c(j,i)/c(i,i)*c(i,:);endenddisp('消元后的矩阵C');cA=c(:,1:m)B=c(:,m+1)x=[];sum=0;x(m)=B(m)/A(m,m);disp('开始回代');for i=m-1:-1:1for j=i+1:msum=sum+A(i,j)*x(j);endx(i)=(B(i)-sum)/A(i,i);sum=0;enddisp('方程的解x1 x2 x3 ...xm')x（2）迭代法clear;clc;% A=[-0.1 -0.2;-0.1 -0.2;-0.2 -0.2];%课本上jacobi迭代法P41 的例子% B=[0.72;0.83;0.84];A=[2/3 1/3;2/3 1/3;1/3 2/3];%第二次试验的方程组B=[39/3;34/3;26/3];c=[];d=[];e=[];x=zeros(1,3);y=zeros(1,3);y(1)=B(1)-A(1,1)*x(2)-A(1,2)*x(3);y(2)=B(2)-A(2,1)*x(1)-A(2,2)*x(3);y(3)=B(3)-A(3,1)*x(1)-A(3,2)*x(2);c=[c;y(1)];d=[d;y(2)];e=[e;y(3)];while i<20 %abs(x(1)-y(1))>0.001&abs(x(2)-y(2))>0.001&abs(x(3)-y(3))>0.001 x=y;%xy(1)=B(1)-A(1,1)*x(2)-A(1,2)*x(3);c=[c;y(1)];y(2)=B(2)-A(2,1)*x(1)-A(2,2)*x(3);d=[d;y(2)];y(3)=B(3)-A(3,1)*x(1)-A(3,2)*x(2);e=[e;y(3)];i=i+1;endplot(c,'-r');hold on;plot(d,'-k');plot(e,'-b');hold offlegend('x1','x2','x3','Location','NorthWest');。
MATLAB上机实验报告：MATLAB求解线性方程组和矩阵的初等计算MATLAB上机实验报告：MATLAB求解线性方程组和矩阵的初等计算计算机语言类课程实验报告课程名称院系学号MATLAB语言及应用电气信息工程学院实验机房计算机公修实验室三班级机器号实验学时2实验成绩专业电气工程及其自动化姓名任课教师实验日期一．实验名称：MATLAB求解线性方程组和矩阵的初等计算二．实验目的和要求1、掌握利用MATLAB程序编辑器编写应用程序的方法；2、掌握MATLAB 求解线性方程组的方法；3、掌握MATLAB进行矩阵的初等计算的方法三．实验内容教材（《MATLAB及其在理工课程中的应用指南，陈怀琛，西安电子科技大学出版社》）P93-1,2四．实验设计方案（实验步骤或开发过程）1、设a矩阵为各个方程的变量的系数，s为各个变量的列向量，b为等式右边的自然数的列向量，利用s=a\\b,即可求出该方程组的解。
2、利用转置和乘法，以及元素乘法分别算出C1,C2,C3，再通过求逆矩阵inv（）这个函数求出C1,C2,C3的逆矩阵。
五．实验中存在问题及解决办法在第一题中当使用右除时，运行出现错误，只有使用左除，这与矩阵除法的定义有关。
六．实验结果1、该方程组的解为：s=[-1.4841-0.68160.5337-1.2429]2、C1=[19-8230d1=[0.00620.0400-0.010612273-0.00460.01690.0030-385429]0.0 1680.02090.0150]C2=[-1516-2436d2=1.0e+015*63-1793-105[-0.9553-0.2391-0.19970.2700226117-600.96670.24200.2021-0.2732194684-10]-0.4473-0.112 0-0.09350.1264-1.1259-0.2818-0.23530.3182]C3=[51624-26d3=不能求C3的逆矩阵-18-12-1572-2-21108-56]七．附录（源程序清单）第一题：求该方程组的解20xx.09.17clcclearalla=[34-7-12;5-742;108-5;-65-210]b=[4;-3;9;-8]s=a\ \b第二题：求出C1,C2,C3,以及它们的逆矩阵d1,d2,d320xx.09.17clcclearallA=[14813;-36-5-9;2-7-12-8]B=[543-2;6-23-8;-13-97]C1=A*B”C2=A”*BC3=A.*Bd1=inv(C1)d2=inv(C2)[m,n]=size(C3);if m==nd3=inv(C3)elsedisp不能求C3的逆矩阵end扩展阅读：袁越强MATLAB上机实验报告一平顶山学院计算机语言类课程实验报告（一）课程名称院系学号实验日期MATLAB语言及应用电气信息工程学院实验机房专业电气工程及其自动化姓名任课教师王凯实验学时23305班级机器号实验成绩二班一．实验名称：MATLAB求解线性方程组和矩阵的初等计算二．实验目的和要求1、掌握利用MATLAB程序编辑器编写应用程序的方法；2、掌握MATLAB求解线性方程组的方法；3、掌握MATLAB进行矩阵的初等计算的方法三．实验内容1、求线性方程组的解3x4y7z12w45x7y4z2w3x8z5w96x5y2z10w8481315432，B6238，求C1=A*B’；C2=A’*B；C3=A.*B，并592、设A36271281397求它们的逆阵。
实验三使用MATLAB解决线性代数问题学院：数计学院班级：1003班姓名：黄晓丹学号：1051020144实验目的：学习MATLAB有关线性代数运算的指令，主要学习运用MATLAB解决矩阵除法，线性方程组的通解，矩阵相似对角化问题，以及解决投入产出分析等应用问题。
实验内容：矩阵转置：A=[1 2;3 4];B=[4 3;2 1];>> A',B'ans =1 32 4ans =4 33 1矩阵加减：A-Bans=-3 -11 3矩阵乘法：A*B,A.*B(数组乘法)||比较矩阵乘法与数组乘法的区别ans=8 520 13ans=4 66 4矩阵除法：A\B,B./Aans=-6 -55 4ans=4 1.50.6667 0.25特殊矩阵生成：zeros（m，n）||生成m行n列的矩阵ones（m，n）||生成m行n列的元素全为一的矩阵eye（n）||生成n阶单位矩阵rand（m，n）||生成m行n列[0 ，1]上均匀分布随机数矩阵zeros(2,3)ans =0 0 00 0 0>> ones(3,3)ans =1 1 11 1 11 1 1>> eye(3)ans =1 0 00 1 00 0 1>> rand(2,4)ans =Columns 1 through 30.9501 0.6068 0.89130.2311 0.4860 0.7621Column 40.45650.0185矩阵处理：trace（A）||返回矩阵的迹diag（A）||返回矩阵对角线元素构成的向量tril（A）||提取矩阵的下三角部分triu（A）||提取矩阵的上三角部分flipud（A）||矩阵上下翻转fliplr（A）||矩阵左右翻转reshape(A,m,n)||将矩阵的元素重排成m行n列矩阵A=[1 2 3;4 5 6;7 8 9];>> t=trace(A),d=diag(A),u=triu(A)t =15d =159u =1 2 30 5 60 0 9 flipud(A),fliplr(A)ans =7 8 94 5 61 2 3 ans =3 2 16 5 49 8 7矩阵特征值与标准型：[V,D]=eig（A）||返回矩阵特征值与特征向量[V J]=Jordan(A)||返回矩阵的相似变换矩阵和若尔当标准型A=[1 2;3 4];>> [V,D]=eig(A)V =-0.8246 -0.41600.5658 -0.9094D =-0.3723 00 5.3723>> [V,J]=jordan(A)V =0.2389 0.76110.5222 -0.5222J =5.3723 00 -0.3723线性方程组求解A=[1 2 1;3 -2 1];B=[1;4];x=A\B x =1.2500 ||求一特解-0.1250>> A=[1 2;3 -2;1 -1];B=[1;4;2];x=A\Bx = ||求得一最小二乘近似解1.2838-0.1757：方阵的相似对角化及应用：A=[1 1/4 0;0 1/2 0;0 1/4 1];[P,T]=eig(A) P =1.0000 0 -0.40820 0 0.81650 1.0000 -0.4082T =1.0000 0 00 1.0000 00 0 0.5000求得三个特征值1,1,0.5，对应特征向量（1,0,0），（0,0,1），（-0.4028,0.8165，-0.4082），由于三个特征向量线性无关，从而A 可相似对角化，即p-1AP=T.那么A∧n=p[1 0 0;0 1 0;0 0 0]p-1，计算的P*diag（[1,1,0]）*inv（P）ans =1.0000 0.50000 00 0 00 0.5000 1.0000所以得到近似解。
实验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 由最简行阶梯型矩阵，得化简后的方程⎪⎩⎪⎨⎧=+--=+--=-+-02200432143214321x x x x x x x x x x x x取x2，x4为自由未知量，扩充方程组为即提取自由未知量系数形成的列向量为基础解系，记所以齐次方程组的通解为解法2： clearA=[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⎩⎨⎧=-=-004321x x x x ⎪⎪⎩⎪⎪⎨⎧====44432221x x x x x x x x ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡+⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡11000011424321x x x x x x ,00111⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε,11002⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε2211εεk k x +=易见，可直接得基础解系所以齐次方程组的通解为二、非齐次线性方程组的求解 Matlab 命令的基本格式：X ＝A\b %系数阵A 满秩时，用左除法求线性方程组AX ＝b 的解注意：A/B 即为AB -1, 而A\B 即为A -1B.C =[A,b];D ＝rref(C) % 求线性方程组AX ＝b 的特解，即D 的最后一列元素【例2】 求下列非齐次线性方程组的解：,00111⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε,11002⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=ε⎪⎪⎪⎩⎪⎪⎪⎨⎧=+=++=++=++=+150650650651655454354332121x x x x x x x x x x x x x 2211εεk k x +=解： clearA=[5 6 0 0 0;1 5 6 0 0;0 1 5 6 0;0 0 1 5 6;0 0 0 1 5]; b=[1;0;0;0;1];format rational %采用有理数近似输出格式,比较format short 看看x=A\b执行后可得所求方程组的解. 作业：【第一题】 求下列非齐次线性方程组的通解.A=[1 2 3 1;1 4 6 2;2 9 8 3;3 7 7 2] B=[3;2;7;12] format rational x=A\B x =⎪⎪⎩⎪⎪⎨⎧=+++=+++=+++=+++1227737389222643324321432143214321x x x x x x x x x x x x x x x x42/31/2684838239393950-7/3【第二题】计算工资问题一个木工，一个电工，一个油漆工，三个人相互同意彼此装修他们自己的房子。
湖南大学电气与信息工程学院 《数值计算》课程 上机实验报告一. 实验目的：了解gauss 消去法和迭代法matlab 算法实现求任意方程组的根。
二. 实验内容：用gauss 消去法和迭代法求解下列线性方程组：263234323923321321321=++=++=++x x x x x x x x x1.求出gauss 消去法的上三角矩阵和方程组的解321,,x x x ，并在命令窗口显示；2.显示迭代法求解过程中所有结果（,,,,,,,,,321131*********NN N x x x x x x x x x ⋯⋯）要求求解精度达到10^-5.三. 算法介绍或方法基础1) 消去法：消元过程：设0)0(11≠a ，令乘数)0(11)0(11/a a m i i -=，做（消去第i 个方程组的i x ）操作1i m ×第1个方程+第i 个方程（i=2，3，.....n ）则第i 个方程变为1)1(2)1(2...i n in i b x a x a =++ 这样消去第2，3，。
，n 个方程的变元i x 后。
原线性方程组变为：⎪⎪⎪⎩⎪⎪⎪⎨⎧=++=++=++)1()1(2)1(2)1(2)1(22)1(22)0(1)0(11)0(11... . .... ...n n nn n n n n n b x a x a b x a x a b x a x a 这样就完成了第1步消元。
回代过程：在最后的一方程中解出n x ，得：)1()1(/--=n nn n n n a b x再将n x 的值代入倒数第二个方程，解出1-n x ，依次往上反推，即可求出方程组的解： 其通项为3,...1-n 2,-n k /)()1(1)1()1(=-=-+=--∑k kk nk j j k kj k kk a x abx高斯赛德尔迭代法：由雅可比迭代公式可知,在迭代的每一步计算过程中是用()k x的全部分量来计算()1+k x的所有分量,显然在计算第i 个分量()1+k ix 时,已经计算出的最新分量()()1111+-+k i k x ,...,x 没有被利用，从直观上看,最新计算出的分量可能比旧的分量要好些.因此，对这些最新计算出来的第1+k 次近似()1+k x 的分量()1+k jx 加以利用,就得到所谓解方程组的高斯—塞德（Gauss-Seidel ）迭代法.把矩阵A 分解成U L D A --= (6)其中()nn a ,...,a ,a diag D 2211=,U ,L --分别为A 的主对角元除外的下三角和上三角部分,于是,方程组(1)便可以写成 ()b Ux x L D +=-即 22f x B x +=其中()()b L D f ,U L D B 1212---=-= (7)以2B 为迭代矩阵构成的迭代法(公式)()()221f x B x k k +=+ (8)称为高斯—塞德尔迭代法(公式),用 量表示的形式为⎩⎨⎧[],...,,k ,n ,,i x a x a b a xi j n i j )k (j ij )k (j ij i ii)k (i21021111111==∑∑--=-=+=++ (9)由此看出,高斯—塞德尔迭代法的一个明显的优点是,在电算时,只需一组存储单元(计算出()1+k ix 后()k ix 不再使用,所以用()1+k ix 冲掉()k ix ,以便存放近似解.四.程序1)消去法：function x=gauss(A,b)n=length(b);A=[A,b];for k=1:(n-1)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);Aendx=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);end2)迭代法：function EX()a=input('请输入系数矩阵a：');b=input('请输入矩阵b:');N=input('请输入最大迭代次数N：');esp=input('请输入近似解的误差限：');if any(diag(a))==0error('系数矩阵错误，迭代终止！')endD=diag(diag(a));X0=zeros(size(b));x1=0;x2=0;x3=0;X1=[x1;x2;x3];h=inv(D)*b;B=inv(D)*(D-a);B1=triu(B);B2=tril(B);k=1;fprintf('高斯-赛德尔迭代法');fprintf('第0次迭代得：')disp(X1');while k<=Nx1=h(1,1)+B1(1,:)*X0;X1=[x1;x2;x3];x2=h(2,1)+B1(2,:)*X0+B2(2,:)*X1;X1=[x1;x2;x3];x3=h(3,1)+B2(3,:)*X1;X1=[x1;x2;x3];if norm(X1-X0,inf)<espfprintf('已满足误差限。
实验一 线性方程组的求解学号： 姓名： 系别专业：实验目的学会用Matlab 软件求解线性方程组.实验内容线性方程组的求解.实验方法和步骤对于线性方程组：AX = b ，其中A 为此线性方程组的系数矩阵，b 为列向量，矩阵 C 为由 A ，b 构成的此方程组的增广矩阵：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=m mnm m n n m n mn m m n n b a a a b a a a b a a a C b b b b x x x X a a a a a a a a a A 212222211112112121212222111211,,,则线性方程组AX=b 的解可分为以下三种个情况：①当m=n ，且rank(A) = rank(C) = n 时，线性方程组AX = b 有唯一解，可通过以下语句给出线性方程组的解：x = inv(A) * b %或者x = A \ b 也可求出方程组的解 x = inv(sym(A)) * b %得出线性方程组的精确解 例1 求解下面的线性方程组.⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡23452314323112344321X解 可验证rank(A) = rank(C) = 4，从而方程组有唯一解，在命令窗口中输入以下语句即可解出此线性方程组：A = [1 2 3 4;4 3 2 1;1 3 2 4;4 1 3 2];b = [5;4;3;2];x = inv(A) * b %利用x = A \ b 所产生的效果是一样的 x =-1.8000 1.8667 3.8667 -2.1333②当rank(A) = rank(C) = r < n 时，线性方程组 AX = b 有无穷多解.此时可以构造出线性方程组的(n - r)个化零向量i x (i = 1，2，…，n - r)，则原线性方程组对应的齐次线性方程组的解^x 可以由i x 的线性组合来表示，即r n x x x x -+++=r -n 2211^ααα其中， αi （i = 1，2，…，n - r)为任意常数.在MATLAB 语言中可以由函数null 直接求出，其调用格式为：Z = null(A) Z = null(A,'r')例2 求下面的线性方程组的解⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡62312244864211224321X解 用下面语句输入矩阵A 与列向量b ，并输入其增广矩阵 C ，从而判定此线性方程组的可解性：A = [1 2 3 4;2 2 1 1;2 4 6 8;4 4 2 2];b = [1;3;2;6];C = [A b]; [rank(A),rank(C)] ans =2 2由于矩阵 A ，C 的秩相同，都等于2，小于未知数的个数4，所以原方程组有无穷多组解，如需求解原线性方程组，可以先求出化零空间Z ，得出满足方程组的一个特解x0：Z = null(A,'r') Z =2.00003.0000 -2.5000 -3.5000 1.0000 0 0 1.0000 x0 = pinv(A) * b x0 =0.9542 0.7328 -0.0763 -0.2977③若 rank(A)≤rank(C)，则方程组 AX = b 无解.这时只能利用Moore-Penrose 广义逆求解出方程的最小二乘解x = pinv(A)*b ，该解不满足原方程组，只是使误差的范数测度||Ax - b||取最小值.例3 讨论下面线性方程组的求解.⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡43212244864211224321X解 输入以下语句，讨论方程组解的情况.A = [1 2 3 4;2 2 1 1;2 4 6 8;4 4 2 2];b = [1;2;3;4];C = [A b]; [rank(A),rank(C)] ans =2 3由于矩阵A 和C 的秩不相等，故原方程组无解，可以使用pinv 函数求解，从而求出原始方程的最小二乘解为：x=pinv(A)*b x =0.5466 0.4550 0.0443 -0.0473由以上语句可知，该解并不满足原线性方程组： norm(A * x - b)ans =0.4472。
matlab实验报告总结1.求一份matlab的试验报告计算方法试验报告3【实验目的】检查各种数值计算方法的长期行为【内容】给定方程组x'(t)=ay(t),y'(t)=bx(t), x(0)=0, y(0)=b的解是x-y 平面上的一个椭圆，利用你已经知道的算法，取足够小的步长，计算上述方程的轨道，看看那种算法能够保持椭圆轨道不变。
（计算的时间步长要足够多）【实验设计】用一下四种方法来计算：1. Euler法2. 梯形法3. 4阶RK法4. 多步法Adams公式【实验过程】1. Euler法具体的代码如下：clear;a=2;b=1;A=[0 a; -b0];U=[];u(:,1)=[0;b];n=1000000;h=6*pi/n;fori=1:n delta(i)=((u(1,i)/a)^2+(u(2,i)/b)^2)^0.5; u(:,i+1)=u(:,i)+h*A*u(:,i);endt=1:n+1;subplot(1, 2,1);plot(1:n,delta);gridon;subplot(1,2,2);plot(u(1,:),u(2,:));gridon;max(abs(delta-ones(1,length(delta)）））；结果如下：2. 梯形法具体的代码如下：clear;a=2;b=1;A=[0 a; -b 0];U=[];u(:,1)=[0;b];n=300;h=6*pi/n;for i=1:n delta(i)=((u(1,i)/a)^2+(u(2,i)/b)^2)^0.5;v1=u(:,i)+h*A*u(:,i);v2=u(:,i)+h*A*(u(:,i)+v1)/2;1u(:,i+1)=u(:,i)+h*A*(u(:,i)+v2)/2;endt=1:n+1;sub plot(1,2,1);plot(1:n,delta);gridon;subplot(1,2,2)；结果如下 3. 4阶RK法clear;a=2;b=1;A=[0 a; -b 0];U=[];u(:,1)=[0;b];n=70;h=6*pi/n;for i=1:n delta(i)=((u(1,i)/a)^2+(u(2,i)/b)^2)^0.5;k1=A*u(:,i); k2=A*(u(:,i)+h/2*k2); k3=A*(u(:,i)+h*k3); k4=A*(u(:,i)+h*k3); u(:,i+1)=u(:,i)+h/6*(k1+2*k2+2*k3+k4);endt=1:n+1 ;subplot(1,2,1);plot(1:n,delta);gridon;subplot(1,2,2)；结果如下：4. 多步法Adams公式clear;a=2;b=1;A=[0 a; -b 0];U=[];u(:,1)=[0;b];n=200;h=6*pi/n;u(:;2)=u(u,1)+h*A*u(:,1);u(:;3)=u(u,2)+h/2*A*(3*u(:,2)-u(:,1));u(:;4)=u(u,3)+h/12*A*(23*u(:,3)-16*u(:,2)+5*u(:, 1)); delta(1)=((u(1,1)/a)^2+(u(2,1)/b^2)^0.5 delta(2)=((u(1,2)/a)^2+(u(2,2)/b^2)^0.5delta(3)=((u(1,3)/a)^2+(u(2,3)/b^2)^0.5for i=4:n delta(i)=((u(1,i)/a)^2+(u(2,i)/b)^2)^0.5;u(:,i+1)=u(:,i)+h/24*A*(55*u(:,i)-59*u(:,i-1)+37 *u(:,i-1)+37*u(:,i-2)-9*u(:,i-3));endt=1:n+1;sub plot(1,2,1);plot(1:n,delta);gridon;subplot(1,2,2)；结果如下：【实验分析】通过这几种方法对比，发现最为稳定的是多步法Adams公式和4阶RK法，其次是梯形法，而欧拉法最为不稳定。
用matlab解线性方程组2008-04-12,17:00一。
高斯消去法1.顺序高斯消去法直接编写命令文件a=[]d=[]'[n,n]=size(a);c=n+1a(:,c)=d;,for,k=1:n-1a(k+1:n,,k:c)=a(k+1:n,,k:c)-(a(k+1:n,k)/,a(k,k))*a(k,,k:c);,,,,,%消去endx=[0,0,0,0]',,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,%回带x(n)=a(n,c)/a(n,n);for,g=n-1:-1:1x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)end2.列主高斯消去法*由于“[r,m]=max(abs(a(k:n,k)))”返回的行是“k:n,k”内的第几行，所以要通过修正来把m,改成真正的行的值。
该程序只是演示程序，真正机器计算不需要算主元素所在列以下各行应为零的值。
直接编写命令文件a=[]d=[],'[n,n]=size(a);c=n+1a(:,c)=d;,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,%(增广)for,k=1:n-1[r,m]=max(abs(a(k:n,k)));,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,%选主m=m+k-1;,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,%(修正操作行的值),,,if(a(m,k)~=0),,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,if(m~=k)a([k,m],:)=a([m,k],:);,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,%换行enda(k+1:n,,k:c)=a(k+1:n,,k:c)-(a(k+1:n,k)/,a(k,k))*a(k,,k:c);,,,,,,%消去endendx=[0,0,0,0]',,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,,%回带x(n)=a(n,c)/a(n,n);for,g=n-1:-1:1x(g)=(a(g,c)-a(g,g+1:n)*x(g+1:n))/a(g,g)end3．分别用顺序高斯消去法和列主高斯消去法解方程组a*x=d，并比较结果a=[0,1,2,3;9,11,23,34;62.5,23.4,15.5,17.2;192.01,124,25.1,59.3],d=[1;1;1;1]顺序高斯消去法：提示“Warning:,Divide,by,zero.”,x,=NaN,NaN,NaN,NaN 列主高斯消去法：x,=-1.2460,2.8796,5.5206,-4.3069由此可见列主高斯消去法可以解决顺序高斯消去法所不能解决的问题。
