数值分析 追赶法 matlab编程
matlab追赶法解常微分方程
研究领域：数学、计算机科学文章标题：深入探讨matlab追赶法解常微分方程在数学和计算机科学领域中，常微分方程是一个重要且广泛应用的课题。
而matlab追赶法作为常微分方程的求解方法，在实际应用中具有重要意义。
本文将以深度和广度兼具的方式，对matlab追赶法解常微分方程这一主题展开全面评估，并撰写一篇有价值的文章，同时结合个人观点和理解，为读者提供深刻的思考。
一、matlab追赶法解常微分方程简介1.1 matlab追赶法基本原理matlab追赶法，又称托马斯算法，是一种用于求解三对角线性方程组的方法。
在常微分方程的数值解法中，常常会遇到需要求解三对角线性方程组的情况，而matlab追赶法正是针对这一问题而提出的高效算法。
1.2 追赶法在常微分方程求解中的应用常微分方程在实际问题中有着广泛的应用，而求解常微分方程的过程中往往需要用到追赶法。
追赶法不仅可以提高计算效率，还可以有效地解决数值稳定性和精度的问题，因此在工程和科学计算中得到了广泛的应用。
二、深入探讨matlab追赶法解常微分方程2.1 算法实现及优化matlab追赶法的实现涉及到矩阵运算、追赶过程和追赶系数的求解等关键步骤。
如何针对不同类型的方程组进行算法优化，是一个需要深入探讨的问题。
通过优化算法，可以提高追赶法的计算效率和数值稳定性，使其在常微分方程求解中发挥更大的作用。
2.2 算法的数值分析通过数值分析，可以更加深入地了解matlab追赶法在解常微分方程过程中的数值特性。
包括收敛性、稳定性、误差分析等方面，这些都是影响算法性能和应用效果的重要因素，需要进行深入的研究和分析。
三、对matlab追赶法解常微分方程的个人观点和理解3.1 算法的优势与局限性matlab追赶法作为一种高效的求解算法，具有较好的稳定性和精度，特别适合于大规模的常微分方程求解。
但在某些特定问题上，追赶法的适用性和效率仍然存在局限性，需要进行合理的选择和应用。
#数值分析Matlab作业
数值分析编程作业2012年12月第二章14.考虑梯形电阻电路的设计，电路如下：电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组：121232343454565676787822/252025202520252025202520250i i V R i i i i i i i i i i i i i i i i i i i i -=-+-=-+-=-+-=-+-=-+-=-+-=-+=这是一个三对角方程组。
设V=220V ，R=27Ω，运用追赶法，求各段电路的电流量。
Matlab 程序如下：function chase () %追赶法求梯形电路中各段的电流量 a=input('请输入下主对角线向量a='); b=input('请输入主对角线向量b='); c=input('请输入上主对角线向量c='); d=input('请输入右端向量d='); n=input('请输入系数矩阵维数n='); u(1)=b(1); for i=2:nl(i)=a(i)/u(i-1); u(i)=b(i)-c(i-1)*l(i); endy(1)=d(1); for i=2:ny(i)=d(i)-l(i)*y(i-1); endx(n)=y(n)/u(n); i=n-1; while i>0x(i)=(y(i)-c(i)*x(i+1))/u(i); i=i-1;end x输入如下: >> chase请输入下主对角线向量a=[0,-2,-2,-2,-2,-2,-2,-2]; 请输入主对角线向量b=[2,5,5,5,5,5,5,5];请输入上主对角线向量c=[-2,-2,-2,-2,-2,-2,-2,0]; 请输入方程组右端向量d=[220/27,0,0,0,0,0,0,0]; 请输入系数矩阵阶数n=8 运行结果如下：x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477第三章14.试分别用(1)Jacobi 迭代法；(2)Gauss-Seidel 迭代法解线性方程组1234510123412191232721735143231211743511512x x x x x ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥=--⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥---⎣⎦⎣⎦⎣⎦ 迭代初始向量(0)(0,0,0,0,0)T x =。
数值分析算法在matlab中的实现
数值分析matlab实现高斯消元法：function[RA,RB,n,X]=gaus(A,b)B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意：因为RA~=RB，所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.')X=zeros(n,1);C=zeros(1,n+1);for p=1:n-1for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1); endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RB<n，所以此方程组有无穷多解.')endendX列主元消去法：function[RA,RB,n,X]=liezhu(A,b)B=[A b];n=length(b);RA=rank(A);RB=rank(B);zhica=RB-RA;if zhica>0,disp('请注意：因为RA~=RB，所以此方程组无解.')returnendif RA==RBif RA==ndisp('请注意：因为RA=RB=n，所以此方程组有唯一解.')X=zeros(n,1);C=zeros(1,n+1);for p=1:n-1[Y,j]=max(abs(B(p:n,p)));C=B(p,:);B(p,:)=B(j+p-1,:);B(j+p-1,:)=C;for k=p+1:nm=B(k,p)/B(p,p);B(k,p:n+1)=B(k,p:n+1)-m*B(p,p:n+1);endendb=B(1:n,n+1);A=B(1:n,1:n);X(n)=b(n)/A(n,n);for q=n-1:-1:1X(q)=(b(q)-sum(A(q,q+1:n)*X(q+1:n)))/A(q,q);endelsedisp('请注意：因为RA=RB<n，所以此方程组有无穷多解.')endendXJacobi迭代法：例1用jacobi迭代法求解代数线性代数方程组,保留四位有效数字(err＝1e-4)其中A=[8-11;2 10-1;11-5];b=[1;4;3]。
追赶法求解三对角线性方程组
追赶法求解三对角线性方程组一 实验目的利用编程方法实现追赶法求解三对角线性方程组。
二 实验内容1、 学习和理解追赶法求解三对角线性方程组的原理及方法；2、 利用MATLAB 编程实现追赶法；3、 举例进行求解，并对结果进行分。
三 实验原理设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵11222=(1)(n 1)()()a c b a c A a n c b n a n ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎣⎦………… 这种方程组称为三对角线性方程组。
显然，A 是上下半宽带都是1的带状矩阵。
设A 的前n-1个顺序主子式都不为零，根据定理2.5的推论，A 有唯一的Crout 分解，并且是保留带宽的。
其中L 是下三角矩阵，U 是单位上三角矩阵。
利用矩阵相乘法，可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⨯⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦……………得到：由上列各式可以得到L 和U 。
引入中间量y ，令y Ux =，则有：已知L 和d ，可求得y 。
则可得到y 的求解表达式：11/12,3,,()(1)*y()=()[()(1)]/y d l i nm i y i li i di y i di m i y i li==-+=--…1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui=*===≤≤+=+++≤≤-=•=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦……………由y Ux =得：111112221u(n 1)(n 1)(n 1)1(n)(n)u x y u x y x y x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦………… 可得到X 的求解表达式：()()1,2,,1()()u()(1)x n y n i n n x i y i i x i ==--=-+…从而得到Ax=d 的解x 。
matlab追赶法解101阶三对角方程组
在探讨MATLAB追赶法解101阶三对角方程组之前，我们首先需要了解什么是追赶法和什么是三对角方程组。
追赶法又称托马斯算法，是一种用于求解带状矩阵（即只有主对角线和两条相邻的对角线上有非零元素的矩阵）的线性方程组的方法。
而三对角矩阵就是只有主对角线和两条相邻的对角线上有非零元素的矩阵。
在实际应用中，求解带状矩阵的线性方程组是非常常见的，特别是在数值计算和科学工程领域。
现在，让我们深入探讨MATLAB追赶法解101阶三对角方程组的方法和具体步骤。
一、MATLAB追赶法解101阶三对角方程组1. 概念介绍101阶三对角方程组是一个非常大的线性方程组，通常使用传统的高斯消元法来求解会耗费大量的时间和计算资源。
而MATLAB追赶法通过利用三对角矩阵的特殊性质，可以有效地简化计算过程，并且节省大量的内存和计算资源。
2. 追赶法步骤（1）将原方程组化为追赶法所需的形式；（2）利用追赶法求解三对角线性方程组。
二、追赶法求解101阶三对角方程组的实现过程1. 将原方程组化为追赶法所需的形式对于101阶三对角方程组，我们首先需要将其化为追赶法所需的形式。
这个过程涉及到选取合适的追赶元和追赶子以及对原方程组的变形，将其化为追赶法能够直接处理的形式。
2. 利用追赶法求解线性方程组一旦将原方程组化为追赶法所需的形式，我们就可以利用追赶法对其进行求解。
追赶法的核心是通过追赶子的迭代计算，逐步求得线性方程组的解。
在MATLAB中，可以使用内置的追赶法求解函数，也可以编写自定义的追赶法算法来实现对101阶三对角方程组的求解。
三、个人观点和理解在实际工程和科学计算中，追赶法是一种非常有效的求解带状矩阵线性方程组的方法。
对于大规模的三对角方程组，特别是高阶的情况，传统的直接求解方法往往会遇到内存和计算资源的限制，而追赶法能够通过精巧的迭代计算，在保证解的精度的显著提高计算效率。
在MATLAB中，通过调用内置的追赶法函数，可以快速地求解大规模的三对角方程组，极大地方便了工程实践中的数值计算工作。
追赶法求解三对角线性方程组
追赶法求解三对角线性方程组一 实验目的利用编程方法实现追赶法求解三对角线性方程组。
二 实验内容1、 学习和理解追赶法求解三对角线性方程组的原理及方法；2、 利用MATLAB 编程实现追赶法；3、 举例进行求解，并对结果进行分。
三 实验原理设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵11222=(1)(n 1)()()a c b a c A a n c b n a n ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎣⎦………… 这种方程组称为三对角线性方程组。
显然，A 是上下半宽带都是1的带状矩阵。
设A 的前n-1个顺序主子式都不为零，根据定理2.5的推论，A 有唯一的Crout 分解，并且是保留带宽的。
其中L 是下三角矩阵，U 是单位上三角矩阵。
利用矩阵相乘法，可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⨯⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦……………得到：由上列各式可以得到L 和U 。
引入中间量y ，令yUx =，则有：已知L 和d ，可求得y 。
则可得到y 的求解表达式：11/12,3,,()(1)*y()=()[()(1)]/y d l i nm i y i li i di y i di m i y i li==-+=--…1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui=*===≤≤+=+++≤≤-=•=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦……………由y Ux =得：111112221u(n 1)(n 1)(n 1)1(n)(n)u x y u x y x y x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦………… 可得到X 的求解表达式：()()1,2,,1()()u()(1)x n y n i n n x i y i i x i ==--=-+… 从而得到Ax=d 的解x 。
数值分析作业-matlab上机作业
数值分析———Matlab上机作业学院：班级：老师：姓名：学号：第二章解线性方程组的直接解法第14题【解】1、编写一个追赶法的函数输入a，b，c，d输出结果x，均为数组形式function x=Zhuiganfa(a,b,c,d)%首先说明：追赶法是适用于三对角矩阵的线性方程组求解的方法，并不适用于其他类型矩阵。
%定义三对角矩阵A的各组成单元。
方程为Ax=d%b为A的对角线元素(1~n)，a为-1对角线元素(2~n)，c为+1对角线元素(1~n-1)。
% A=[2 -1 0 0% -1 3 -2 0% 0 -2 4 -3% 0 0 -3 5]% a=[-1 -2 -3];c=[-1 -2 -3];b=[2 3 4 5];d=[6 1 -2 1];n=length(b);u(1)=b(1);y(1)=d(1);for i=2:nl(i)=a(i-1)/u(i-1);%先求l(i)u(i)=b(i)-c(i-1)*l(i);%再求u(i)%A=LU,Ax=LUx=d,y=Ux,%Ly=d,由于L是下三角矩阵，对角线均为1，所以可求y(i)y(i)=d(i)-l(i)*y(i-1);endx(n)=y(n)/u(n);for i=(n-1):-1:1%Ux=y,由于U是上三角矩阵，所以可求x(i)x(i)=(y(i)-c(i)*x(i+1))/u(i);end2、输入已知参数>>a=[2 2 2 2 2 2 2];>>b=[2 5 5 5 5 5 5 5];>>c=[2 2 2 2 2 2 2];>>d=[220/27 0 0 0 0 0 0 0];3、按定义格式调用函数>>x=zhuiganfa(a,b,c,d)4、输出结果x=[8.147775166909105 -4.073701092835030 2.036477565178471 -1.017492820111148 0.507254485099400 -0.250643392637350 0.119353996493976 -0.047741598597591]第15题【解】1、编写一个程序生成题目条件生成线性方程组A x=b 的系数矩阵A 和右端项量b ，分别定义矩阵A 、B 、a 、b 分别表示系数矩阵，其中1(10.1;,1,2,...,)j ij i i a x x i i j n -==+=或1(,1,2,...,)1ij a i j n i j ==+-分别构成A 、B 对应右端项量分别a 、b 。
数值分析的MATLAB程序
列主元法function lianzhuyuan(A,b)n=input('请输入n:') %选择阶数A=zeros(n,n); %系数矩阵Ab=zeros(n,1); %矩阵bX=zeros(n,1); %解Xfor i=1:nfor j=1:nA(i,j)=(1/(i+j-1)); %生成hilbert矩阵Aendb(i,1)=sum(A(i,:)); %生成矩阵bendfor i=1:n-1j=i;top=max(abs(A(i:n,j))); %列主元k=j;while abs(A(k,j))~=top %列主元所在行k=k+1;endfor z=1:n %交换主元所在行a1=A(i,z);A(i,z)=A(k,z);A(k,z)=a1;enda2=b(i,1);b(i,1)=b(k,1);b(k,1)=a2;for s=i+1:n %消去算法开始m=A(s,j)/A(i,j); %化简为上三角矩阵A(s,j)=0;for p=i+1:nA(s,p)=A(s,p)-m*A(i,p);endb(s,1)=b(s,1)-m*b(i,1);endendX(n,1)=b(n,1)/A(n,n); %回代开始for i=n-1:-1:1s=0; %初始化sfor j=i+1:ns=s+A(i,j)*X(j,1);endX(i,1)=(b(i,1)-s)/A(i,i);endX欧拉法clcclear% 欧拉法p=10; %贝塔的取值T=10; %t取值的上限y1=1; %y1的初值r1=1; %y2的初值%输入步长h的值h=input('欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代过程for t=0:h:TY=(exp(-t));R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t);y=y1+h*(-y1);y1=y;r=r1+h*(y1-p*r1);r1=r;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=[0:h:T];% 红线为解析解，'x'为数值解plot(t,S1,'r',t,S3,'x')改进欧拉法clcclearp=10; %贝塔的取值T=10; %t取值的上限y1=1; %y1的初值r1=1; %y2的初值%输入步长h的值h=input('改进欧拉法please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代过程for t=0:h:TY=(exp(-t));R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t);k1=-y1;l1=y1-p*r1;k2=-(y1+h*k1);l2=y1+h*k1-p*(r1+h*l1);y=y1+h*(k1+k2)/2;r=r1+h*(k1+k2)/2;r1=r;y1=y;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=[0:h:T];% 红线为解析解，'x'为数值解plot(t,S1,'r',t,S3,'x')高斯-赛德尔function gaosisaideern=input('n='); %阶数tol=input('tol='); %迭代精度A=zeros(n,n);b=zeros(n,1); %生成b向量for i=1:n %给Hilbert矩阵和b向量赋值for j=1:nA(i,j)=(1/(i+j-1));endb(i,1)=sum(A(i,:));endy=zeros(n,1); %迭代解x1=zeros(n,1); %准确解for i=1:ny(i,1)=0; %迭代解赋初值x1(i,1)=1; %生成准确解endk=0;while norm(y-x1)>=tol %精度控制(采用自动步数控制) k=k+1;for i=1:n %迭代开始a1=0;a2=0;for j=1:i-1a1=a1+A(i,j)*y(j,1);endfor j=i+1:na2=a2+A(i,j)*y(j,1);endy(i,1)=(b(i,1)-a1-a2)/A(i,i);endenddisp('迭代步数k')kdisp(y) %显示yend最速下降法function gaosisaideern=input('阶数n='); %阶数tol=input('迭代精度tol='); %迭代精度eps=input('最速下降法eps=');A=zeros(n,n);b=zeros(n,1); %生成b向量for i=1:n %给Hilbert矩阵和b向量赋值for j=1:nA(i,j)=(1/(i+j-1));endb(i,1)=sum(A(i,:));endy=zeros(n,1); %迭代解x1=zeros(n,1); %准确解t=zeros(n,1);r=zeros(n,1);for i=1:ny(i,1)=0; %迭代解赋初值x1(i,1)=1; %生成准确解endr=b-A*y;while norm(r)>=eps; %先进行最速下降法求得进行赛德尔迭代的初始解yt=(r'*r)/(r'*A*r);s1=t*r;y=y+s1;r=b-A*y;endk=0;while norm(y-x1)>=tol %精度控制(采用自动步数控制)k=k+1;for i=1:n %迭代开始a1=0;a2=0;for j=1:i-1a1=a1+A(i,j)*y(j,1);endfor j=i+1:na2=a2+A(i,j)*y(j,1);endy(i,1)=(b(i,1)-a1-a2)/A(i,i);endenddisp('迭代步数k')disp(k)disp(y) %显示y四阶龙格-库塔法clcclearp=10; %贝塔的取值T=10; %t取值的上限y1=1; %y1的初值r1=1; %y2的初值%输入步长h的值h=input('四阶龙格please input number(h=1 0.5 0.25 0.125 0.0625):h=') ; if h>1 or h<0breakendS1=0:T/h;S2=0:T/h;S3=0:T/h;S4=0:T/h;i=1;% 迭代过程for t=0:h:TY=(exp(-t));R=(1/(p-1))*exp(-t)+((p-2)/(p-1))*exp(-p*t);k1=-y1;l1=y1-p*r1;k2=-(y1+h*k1/2);l2=y1+h*k1/2-p*(r1+h*l1/2);k3=-(y1+h*k2/2);l3=y1+h*k2/2-p*(r1+h*l2/2);k4=-(y1+h*k3);l4=y1+h*k3-p*(r1+h*l3);y=y1+h*(k1+2*k2+2*k3+k4)/6;r=r1+h*(l1+2*l2+2*l3+l4)/6;r1=r;y1=y;S1(i)=Y;S2(i)=R;S3(i)=y;S4(i)=r;i=i+1;endt=[0:h:T];% 红线为解析解，'x'为数值解plot(t,S1,'r',t,S3,'x')。
追赶法求解三对角线性方程组
追赶法求解三对角线性方程组一 实验目的利用编程方法实现追赶法求解三对角线性方程组。
二 实验内容1、 学习和理解追赶法求解三对角线性方程组的原理及方法；2、 利用MA TLAB 编程实现追赶法；3、 举例进行求解，并对结果进行分。
三 实验原理设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵11222=(1)(n 1)()()a c b a c A a n c b n a n ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎣⎦………… 这种方程组称为三对角线性方程组。
显然，A 是上下半宽带都是1的带状矩阵。
设A 的前n-1个顺序主子式都不为零，根据定理2.5的推论，A 有唯一的Crout 分解，并且是保留带宽的。
其中L 是下三角矩阵，U 是单位上三角矩阵。
利用矩阵相乘法，可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⨯⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦……………得到：由上列各式可以得到L 和U 。
引入中间量y ，令yUx =，则有：已知L 和d ，可求得y 。
则可得到y 的求解表达式：11/12,3,,()(1)*y()=()[()(1)]/y d l i nm i y i li i di y i di m i y i li==-+=--…1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui=*===≤≤+=+++≤≤-=∙=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦……………由y Ux =得：111112221u(n 1)(n 1)(n 1)1(n)(n)u x y u x y x y x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦………… 可得到X 的求解表达式：()()1,2,,1()()u()(1)x n y n i n n x i y i i x i ==--=-+… 从而得到Ax=d 的解x 。
数值分析列主元高斯消去顺序高斯平方根法追赶法
课题名称：课题一解线性方程组的直接方法解决的问题：给定三个不同类型的线性方程组，用适当的直接法求解。
采用的数值方法：对第一个普通的线性方程组，采用了高斯顺序消去法和高斯列主元消去法。
对第二个正定线性方程组，采用了平方根法。
对第三个三对角线性方程组，采用了追赶法。
算法程序：(1)普通的线性方程组①顺序消去法#include<stdio.h>#include<math.h>int main(void){float A[10][10]= {{4,2,-3,-1,2,1,0,0,0,0},{8,6,-5,-3,6,5,0,1,0,0},{4,2,-2,-1,3,2,-1,0,3,1},{0,-2,1,5,-1,3,-1,1,9,4},{-4,2,6,-1,6,7,-3,3,2,3},{8,6,-8,5,7,17,2,6,-3,5},{0,2,-1,3,-4,2,5,3,0,1},{16,10,-11,-9,17,34,2,-1,2,2},{4,6,2,-7,13,9,2,0,12,4},{0,0,-1,8,-3,-24,-8,6,3,-1}};float b[10]= {5,12,3,2,3,46,13,38,19,-21}; float x[10]= {0};float Aik,S,temp;int i,j,k;int size=10;for(k=0; k<size-1; k++){if(!A[k][k])return -1;for(i=k+1; i<size; i++){Aik=A[i][k]/A[k][k];for(j=k; j<size; j++){A[i][j]=A[i][j]-Aik*A[k][j]; }b[i]=b[i]-Aik*b[k];}}printf("A[]\n");for(i=0; i<size; i++){for(j=0; j<size; j++)printf("%f ",A[i][j]);printf("\n");}printf("b[]\n");for(i=0; i<size; i++)printf("%f ",b[i]);printf("\n\n");x[size-1]=b[size-1]/A[size-1][size-1]; for(k=size-2; k>=0; k--){S=b[k];for(j=k+1; j<size; j++){S=S-A[k][j]*x[j];}x[k]=S/A[k][k];}printf("x[]=\n");for(i=0; i<size; i++)printf("%f ",x[i]);return 0;}②列主元消去法#include<stdio.h>#include<math.h>int main(void){float A[10][10]= {{4,2,-3,-1,2,1,0,0,0,0}, {8,6,-5,-3,6,5,0,1,0,0},{4,2,-2,-1,3,2,-1,0,3,1},{0,-2,1,5,-1,3,-1,1,9,4},{-4,2,6,-1,6,7,-3,3,2,3},{8,6,-8,5,7,17,2,6,-3,5},{0,2,-1,3,-4,2,5,3,0,1},{16,10,-11,-9,17,34,2,-1,2,2},{4,6,2,-7,13,9,2,0,12,4},{0,0,-1,8,-3,-24,-8,6,3,-1}};float b[10]= {5,12,3,2,3,46,13,38,19,-21}; float x[10]= {0};float Aik,S,temp;int i,j,k;float max;int col;int size=10;for(k=0; k<size-1; k++){max=fabs(A[k][k]);col=k;for(i=k; i<size; i++){if(max<fabs(A[i][k])) {max=fabs(A[i][k]); col=i;}}for(j=k; j<size; j++){temp=A[col][j];A[col][j]=A[k][j];A[k][j]=temp;}temp=b[col];b[col]=b[k];b[k]=temp;if(!A[k][k])return -1;for(i=k+1; i<size; i++){Aik=A[i][k]/A[k][k]; for(j=k; j<size; j++){A[i][j]=A[i][j]-Aik*A[k][j]; }b[i]=b[i]-Aik*b[k];}}printf("A[]\n");for(i=0; i<size; i++){for(j=0; j<size; j++)printf("%f ",A[i][j]);printf("\n");}printf("b[]\n");for(i=0; i<size; i++)printf("%f ",b[i]);printf("\n\n");x[size-1]=b[size-1]/A[size-1][size-1]; for(k=size-2; k>=0; k--){S=b[k];for(j=k+1; j<size; j++){S=S-A[k][j]*x[j]; }x[k]=S/A[k][k];}printf("x[]=\n");for(i=0; i<size; i++)printf("%f ",x[i]);return 0;}(2)对称正定线性方程组平方根法：#include <stdio.h>#include <math.h>#define n 8int main(void){float A[8][8]={{4,2,-4,0,2,4,0,0},{2,2,-1,-2,1,3,2,0},{-4,-1,14,1,-8,-3,5,6},{0,-2,1,6,-1,-4,-3,3},{2,1,-8,-1,22,4,-10,-3},{4,3,-3,-4,4,11,1,-4},{0,2,5,-3,-10,1,14,2},{0,0,6,3,-3,-4,2,19}};float g[8][8]= {0};float b[8]= {0,-6,6,23,11,-22,-15,45}; float x[8]= {0};float y[8]= {0};int k,m,i,sq;for(k=0; k<n; k++){float p=0,q=0,s=0;for(m=0; m<=k-1; m++){p=p+A[k][m]*A[k][m];}g[k][k]=sqrt(A[k][k]-p);A[k][k]=g[k][k];for(i=k+1; i<n; i++){q=0;for(m=0; m<=k-1; m++){q=q+A[i][m]*A[k][m];}g[i][k]=(A[i][k]-q)/A[k][k]; A[i][k]=g[i][k];}s=0;for(m=0; m<=k-1; m++){s=s+A[k][m]*y[m];}y[k]=(b[k]-s)/A[k][k];}x[n-1]=y[n-1]/A[n-1][n-1];for(k=n-2; k>=0; k--){float sum=0;for(m=k+1; m<n; m++){sum=sum+A[m][k]*x[m];}x[k]=(y[k]-sum)/A[k][k];}for(sq=0; sq<n; sq++){printf("%f ",x[sq]);}return 0;}(3)三对角线性方程组追赶法#include <stdio.h>#include <math.h>#define n 10int main(void){float a[10]={4,4,4,4,4,4,4,4,4,4};float c[9]={-1,-1,-1,-1,-1,-1,-1,-1,-1};float d[9]={-1,-1,-1,-1,-1,-1,-1,-1,-1}; float b[10]={7,5,-13,2,6,-12,14,-4,5,-5}; float x[10]={0};float y[10]={0};float arf[10]={0};float bt[9]={0};arf[0]=a[0];int i;for(i=0;i<n-1;i++){bt[i]=c[i]/arf[i];arf[i+1]=a[i+1]-d[i+1]*bt[i];//printf("%f %f \n",bt[i],arf[i+1]); }y[0]=b[0]/arf[0];//printf("%f\n",y[0]);for(i=1;i<n;i++){y[i]=(b[i]-d[i]*y[i-1])/arf[i];}//printf("%f\n",y[1]);x[n-1]=y[n-1];for(i=n-2;i>=0;i--){x[i]=y[i]-bt[i]*x[i+1]; }for(i=0;i<n;i++)printf("%lf ",x[i]);return 0;}数值结果：(1) 普通的线性方程组①顺序消去法②列主元消去法(2) 对称正定线性方程组平方根法：(3)三对角线性方程组追赶法：对实验计算结果的讨论和分析：(1) 普通的线性方程组①顺序消去法x1~x10的绝对误差：0.000001，-0.000001,0.000001,0,0.000001,0,0.000002,0,0,0x1~x10的相对误差：0.000001，0.000001,-1,0,0.0000005,0,0.00000067,0,0,0误差很小，基本可以忽略。
数值分析实验报告三
grid
[k,x,wuca,yx]=erfen (﹣1,1,10^-5)
2）运行结果
ans =
0 -1.0000 1.0000 0 1.0000 -11.6321 10.7183 -1.0000
ans =
1.0000 0 1.0000 0.5000 0.5000 -1.0000 10.7183 4.6487
ans =
11.0000 0.0898 0.0908 0.0903 0.0005 -0.0076 0.0033 -0.0021
ans =
12.0000 0.0903 0.0908 0.0906 0.0002 -0.0021 0.0033 0.0006
ans =
13.0000 0.0903 0.0906 0.0905 0.0001 -0.0021 0.0006 -0.0008
ans =
7.0000 0.1256 0.0008 0.0033 0.0262
ans =
8.0000 0.1240 0.0002 0.0016 0.0129
ans =
9.0000 0.1233 0.0000 0.0007 0.0056
ans =
9.0000 0.1233 0.0000 0.0007 0.0056
(2)、Use the iteration method ,the initial value .
2、The equation has two roots near 0.1.
Determine them by means ofNewton’s method.
(with accuracy )
3、用迭代法求方程 附近的一个根。方程写成下
k = 9
数值分析MATLAB科学计算—线性方程组
科学计算—理论、方法及其基于ＭＡＴＬＡＢ的程序实现与分析 三、 解线性方程组（线性矩阵方程）解线性方程组是科学计算中最常见的问题。
所说的“最常见”有两方面的含义：１） 问题的本身是求解线性方程组；２） 许多问题的求解需要或归结为线性方程组的求解。
关于线性方程组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注：数值稳定的算法高斯列主元素消去法就是在消元的每一步选取（列）主元素—一列中绝对值最大的元取做主元素，高斯列主元素消去法是数值稳定的方法。
数值分析编程总结
1。
LU分解:[L U]=lu（A）;2. 追赶法function x=zhuiganfa（A,b）[n,n]=size（A）;for i=1：nif(i==1）l(i)=A(i，i）;y(i）=b(i)/l(i）；elsel(i)=A(i,i）—A（i，i—1)＊u(i—1）;y(i）=（b（i）—y（i—1）*A（i，i—1))/l(i）；endif（i〈n)u(i)=A（i，i+1）/l（i)；endendx（n）=y(n）for j=n—1:-1:1x(j)=y（j）—u（j）*x（j+1）；end数值试验：n=101;a= 12 1 0……………………………….。
…。
01 12 1 00 1 12 1 ...........................。
. 00 0 1 12 1 011 12关键是如何定义上述矩阵：>〉n=101；c1=ones(1，n-1）;a1=diag(c1，—1)；这个—1说明行位置—1c2=12＊ones(1,n);a2=diag(c2）；c3=ones(1,n—1）;a3=diag（c3,1)；a=a1+a2+a3；3. 拉格朗日插值function yh=lage（x,y，xh)n=length(x);m=length（xh);yh=zeros(1，m)；c1=ones（n—1,1）；c2=ones（1,m）；for i=1：nxp=x（[1:i—1 i+1：n]）;yh=yh+y(i)＊prod（（c1＊xh-xp’*c2)。
/（x(i)-xp’*c2)）；endend>> x=［11，12］；〉> y=[2，4]；〉〉xh=［11.75];>> lage（x,y,xh)ans =3。
50004 最小二乘法1.2.编程函数为：function z = erchen(x,y）x1=ones（5,1）；A=[x1，x,-x.*y]; 注意点乘z=A\y；注意左除a=z（1）;b=z(2)；c=z(3);end输入:≻≻ x=[1953 1964 1982 1990 2000］’；≻≻ y=［5。
《数值分析》上机实验报告
数值分析上机实验报告《数值分析》上机实验报告1.用Newton 法求方程 X 7-X 4+14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
1.1 理论依据：设函数在有限区间[a ，b]上二阶导数存在，且满足条件{}αϕ上的惟一解在区间平方收敛于方程所生的迭代序列迭代过程由则对任意初始近似值达到的一个中使是其中上不变号在区间],[0)(3,2,1,0,)(')()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20)()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f ab c f x f b a x f b f x f k k k k k k ==-==∈≤-≠>+令)9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3225333647>⋅''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f故以1.9为起点⎪⎩⎪⎨⎧='-=+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。
当前后两个的差<=ε时，就认为求出了近似的根。
本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。
1.2 C语言程序原代码：#include<stdio.h>#include<math.h>main(){double x2,f,f1;double x1=1.9; //取初值为1.9do{x2=x1;f=pow(x2,7)-28*pow(x2,4)+14;f1=7*pow(x2,6)-4*28*pow(x2,3);x1=x2-f/f1;}while(fabs(x1-x2)>=0.00001||x1<0.1); //限制循环次数printf("计算结果：x=%f\n",x1);}1.3 运行结果：1.4 MATLAB上机程序function y=Newton(f,df,x0,eps,M)d=0;for k=1:Mif feval(df,x0)==0d=2;breakelsex1=x0-feval(f,x0)/feval(df,x0);ende=abs(x1-x0);x0=x1;if e<=eps&&abs(feval(f,x1))<=epsd=1;breakendendif d==1y=x1;elseif d==0y='迭代M次失败';elsey= '奇异'endfunction y=df(x)y=7*x^6-28*4*x^3;Endfunction y=f(x)y=x^7-28*x^4+14;End>> x0=1.9;>> eps=0.00001;>> M=100;>> x=Newton('f','df',x0,eps,M);>> vpa(x,7)1.5 问题讨论：1.使用此方法求方解，用误差来控制循环迭代次数，可以在误差允许的范围内得到比较理想的计算结果。
数值分析中求解线性方程组的MATLAB程序(6种)
数值分析中求解线性方程组的MATLAB程序（6种）1.回溯法（系数矩阵为上三角）function X=uptrbk(A,B)%求解方程组，首先化为上三角，再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));2.系数矩阵为下三角function x=matrix_down(A,b)%求解系数矩阵是下三角的方程组n=length(b);x=zeros(n,1);x(1)=b(1)/A(1,1);for k=2:1:nx(k)=(b(k)-A(k,1:k-1)*x(1:k-1))/A(k,k);end3.普通系数矩阵（先化为上三角，在用回溯法）function X=uptrbk(A,B)%求解方程组，首先化为上三角，再调用函数求解[N,N]=size(A);X=zeros(N,1);C=zeros(1,N+1);Aug=[A B];for p=1:N-1[Y,j]=max(abs(Aug(p:N,p)));C=Aug(p,:);Aug(p,:)=Aug(j+p-1,:);Aug(j+p-1,:)=C;if Aug(p,p)==0'A was singular.No unique solution.'break;endfor k=p+1:Nm=Aug(k,p)/Aug(p,p);Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);endendD=Aug;X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));4.三角分解法function [X,L,U]=LU_matrix(A,B)%A是非奇异矩阵%AX=B化为LUX=B，L为下三角，U为上三角%程序中并没有真正解出L和U，全部存放在A中[N,N]=size(A);X=zeros(N,1);Y=zeros(N,1);C=zeros(1,N);R=1:N;for p=1:N-1[max1,j]=max(abs(A(p:N,p)));C=A(p,:);A(p,:)=A(j+p-1,:);A(j+p-1,:)=C;d=R(p);R(p)=R(j+p-1);R(j+p-1)=d;if A(p,p)==0'A is singular.No unique solution'break;endfor k=p+1:Nmult=A(k,p)/A(p,p);A(k,p)=mult;A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);endendY(1)=B(R(1));for k=2:NY(k)=B(R(k))-A(k,1:k-1)*Y(1:k-1);endX(N)=Y(N)/A(N,N);for k=N-1:-1:1X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);endL=tril(A,-1)+eye(N)U=triu(A)5.雅克比迭代法function X=jacobi(A,B,P,delta,max1);%雅克比迭代求解方程组N=length(B);for k=1:max1for j=1:NX(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);enderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)breakendendX=X';k6.盖斯迭代法function X=gseid(A,B,P,delta,max1);%盖斯算法，求解赋初值的微分方程N=length(B);for k=1:max1for j=1:Nif j==1X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1);elseif j==NX(N)=(B(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);elseX(j)=(B(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);endenderr=abs(norm(X'-P));relerr=err/(norm(X)+eps);P=X';if (err<delta)|(relerr<delta)break;endendX=X';k。
数值分析中常用的matlab程序
1.% 最小二乘法拟合数据点方法1：% 左除右除：xA=B ==> x=B/A | Ax=B ==> x=B\A% A 表示拟合函数的组合，如：多项式插值，A=[1,x,x.^2,...,x.^n],表示拟合函数为% 多项式：s(x)=a0+a1*x+...+an*x^n;又如：A=[log(x),cos(x),exp(x)]%则表示拟合函数为s(x)=a0*ln(x)+a1*cos(x)+a2*exp(x)% 法方程为：A'*A*z=A'*y ==> A*z=y ==> z=A\y z=(a0,...,an)'% Date: 2012-1-1clear;clc;x=[0 0.25 0.5 0.75 1]';y=[1 1.284 1.6487 2.1170 2.7183]';%索要拟合的数据点x1=ones(size(x),1);A=[x1 x x.^2];%拟合函数Z=A\y %A中每列函数的参数% 最小二乘法拟合数据方法2：采用polyfit函数% Date:2012-1-1clear;clc;x=[0 0.25 0.5 0.75 1]';y=[1 1.284 1.6487 2.1170 2.7183]';%索要拟合的数据点p=polyfit(x,y,2) %polyfit(x,y,n),(x,y)为数据点坐标，n为拟合多项式阶数，p 为% 所求拟合多项式的幂次从高到低排列的系数。
2.% 复合梯形公式计算积分值% 输入：fun--积分函数；a，b--积分区间；n--区间等分数% 输出：I--数值积分结果% 调用格式(ex)：re=ftrapz(@fun1,0,1,10)% 2012-1-1function I=ftrapz(fun,a,b,n)h=(b-a)/n;%区间等分x=linspace(a,b,n+1);%将a到b的区间等分成(n+1)-1个区间，数据点有(n+!)个y=feval(fun,x);I=h*(0.5*y(1)+sum(y(2:n))+0.5*y(n+1));%积分原函数% Date：2012-1-1function y=fun1(x)y=exp(-x);3.% 复合simpson公式求积分% 输入：fun--积分函数；a，b--积分区间；n--区间等分数% 输出：I--数值积分结果% 调用格式(ex)：re=fsimpson(@fun1,0,1,10)% 2012-1-1function I=fsimpson(fun,a,b,n)h=(b-a)/n;x=linspace(a,b,2*n+1);y=feval(fun,x);I=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1));4.% 两点GS-Legendre公式求积分% 输入：fun--积分函数；a，b--积分区间；% 输出：I--数值积分结果% 调用格式(ex)：re=GSLege(@fun2,0,1)% 2012-1-1function I=GSLege(fun,a,b)%将区间[a.b]通过变量替换x=(a+b)/2-(b-a)/2*t变到[-1,1]%其中t取GS-Lege的高斯点m1=feval(fun,(a+b)/2+(b-a)/2*(-1/sqrt(3)));m2=feval(fun,(a+b)/2+(b-a)/2*(1/sqrt(3)));I=(b-a)/2*(m1+m2);% GS-Legendre公式的积分函数% Date:2012-1-1function y=fun2(x)y=sin(x)/x;5.% 追赶法求解线性方程组Ax=b，其中A是三对角方阵%function x=tridiagsolver(A,b)clear;clc;A=[2 -1 0 0;-1 3 -2 0;0 -2 4 -3;0 0 -3 5];%三对角矩阵，线性方程组系数矩阵b=[6,1,-2,1]';%[n,n]=size(A);for i=1:nif(i<2)l(i)=A(i,i);y(i)=b(i)/l(i);u(i)=A(i,i+1)/l(i);elseif i<nl(i)=A(i,i)-A(i,i-1)*u(i-1);y(i)=(b(i)-y(i-1)*A(i,i-1))/l(i);u(i)=A(i,i+1)/l(i);elsel(i)=A(i,i)-A(i,i-1)*u(i-1);y(i)=(b(i)-y(i-1)*A(i,i-1))/l(i);endendx(n)=y(n);for j=n-1:-1:1x(j)=y(j)-u(j)*x(j+1);endx6.% SOR迭代求解非线性方程组Ax=b% 输入：A--系数矩阵；b--；omega--松弛因子（0~2);tol--精度% 输出：x--方程的解向量;iter--迭代次数;% 调用格式(ex)：[x,iter=sor(A,b,1.1,1e-4)% 2012-1-1%function [x,iter]=sor(A,b,omega,tol)%{%}clear;clc;A=[2 -1 0;-1 3 -1;0 -1 2];b=[1 8 -5]';omega=1.1;tol=1e-4;D=diag(diag(A));%diag(A)返回的是A的对角元组成的列向量；diag(b)返回的是以列向量b为对角元的方阵；L=D-tril(A);%tril(A)返回的是矩阵A的下三角矩阵；U=D-triu(A);%triu(A)返回的是矩阵A的上三角矩阵；x=zeros(size(b));%给定迭代初始值零向量iter=1;while iter<500x=(D-omega*L)\(omega*b+(1-omega)*D*x+omega*U*x);%迭代格式error=norm(b-A*x)/norm(b);%收敛条件if error<tolbreak;enditer=iter+1;endif iter>= 500fprintf('root not found!');end7.% 牛顿法求解非线性方程的根% 输入：fun--非线性函数；dfun--非线性函数导数；x0--初始值；tol--精度；% 输出：x--非线性方程数值根% 调用格式(ex)：x=newton(@fun3,@dfun3,6,1e-3)% 2012-1-1function x=newton(fun,dfun,x0,tol)iter=1;if abs(feval(fun,x0))<tolx=x0;return;endwhile iter<500if iter==1x=x0-feval(fun,x0)/feval(dfun,x0);if abs(feval(fun,x))<tolbreak;endelsex=x-feval(fun,x)/feval(dfun,x);if abs(feval(fun,x))<tolbreak;endenditer=iter+1;endif iter==500fprintf('not successful!');x=NaN;end% newton的函数文件% Date：2012-1-1function y=fun3(x)y=3.*x.^3-8.*x.^2-8.*x-11;% newton的导函数文件% Date：2012-1-1function y=dfun3(x)y=9.*x.^2-16.*x-8;8.% 两点割线法求解非线性方程的根% 输入：fun--非线性函数；a,b--两个初始值；tol--精度；% 输出：x--非线性方程数值根% 调用格式(ex)：x=gexian(@fun3,3,6,1e-3)% 2012-1-1function x=gexian(fun,a,b,tol)iter=1;xk1=a;xk2=b;while iter<500xk3=xk2-feval(fun,xk2)*(xk2-xk1)/(feval(fun,xk2)-feval(fun,xk1));if abs(feval(fun,xk3))<tolbreak;elsexk1=xk2;xk2=xk3;enditer=iter+1;endif iter==500fprintf('not successful!');x=NaN;elsex=xk3;end9.% 乘幂法求矩阵的按模最大特征值及其特征向量% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；% 输出：x--特征向量；lam--按模最大特征值% 调用格式(ex)：[lam,x]=power(A,v0,1e-2)% 2012-1-1function [lam,x]=matrixpower(A,v0,tol)v1=A*v0;v2=A*v1;sum=0;p=0;for i=1:size(v1)if v1(i)~=0sum=sum+v2(i)/v1(i);p=p+1;endendlam0=sum/p;iter=2;vk1=v2;while iter<500vk2=A*vk1;sum=0;p=0;for i=1:size(vk1)if vk1(i)~=0sum=sum+vk2(i)/vk1(i);p=p+1;endendlam=sum/p;if abs(lam-lam0)<tolbreak;elselam0=lam;vk1=vk2;endendif iter==500fprintf('not successful!');lam=NaN;elsex=vk2;end9_2% 改进乘幂法求矩阵的按模最大特征值及其特征向量% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；% 输出：x--特征向量；lam--按模最大特征值% 调用格式(ex)：[lam,x]=pMatrixPower(A,v0,1e-2)% Date:2012-1-2function [lam,x]=pMatrixPower(A,v0,tol)[ty,ti]=max(abs(v0));%返回v0中元素绝对值最大的元素值与下标,ti为下标lam0=v0(ti);u0=v0/lam0;iter=1;while iter<500v1=A*u0;[tv,ti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;if abs(lam0-lam1)<tolbreak;elselam0=lam1;iter=iter+1;endendif iter>=500fprintf('not successful!');lam=NaN;elselam=lam1;x=u0;end10.% 反幂法求矩阵的按模最小特征值及其特征向量% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；% 输出：x--特征向量；lam--按模最大特征值% 调用格式(ex)：[lam,x]=invMatrixPower(A,v0,1e-2)% Date:2012-1-2function [lam,x]=invMatrixPower(A,v0,tol)[ty,ti]=max(abs(v0));lam0=v0(ti);u0=v0/lam0;iter=1;while iter<500v1=A\u0;[tv,ti]=max(abs(v1));lam1=v1(ti);u0=v1/lam1;if abs(1/lam0-1/lam1)<tolbreak;elselam0=lam1;iter=iter+1;endendif iter>=500fprintf('not successful!');lam=NaN;elselam=1/lam1;x=u0;end11.% 欧拉方法求解一阶常微分方程初值问题% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;% 输出：y--常微分方程在区间[a,b]上各点的数值解；% 调用格式(ex)：y=odeEuler(@fun4,0,1,1,10)% Date:2012-1-2function y=odeEuler(fun,a,b,y0,n)h=(b-a)/n;y(1)=y0;x=a:h:b;for i=1:ny(i+1)=y(i)+h*feval(fun,x(i),y(i));end% 常微分方程fun4=0% Date：2012-1-2function re=fun4(x,y)re=y-2*x/y;12.% 改进欧拉公式（预估--校正）方法求解一阶常微分方程初值问题% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;% 输出：y--常微分方程在区间[a,b]上各点的数值解；% 调用格式(ex)：y=pOdeEuler(@fun4,0,1,1,10)% Date:2012-1-2function y=pOdeEuler(fun,a,b,y0,n)h=(b-a)/n;y(1)=y0;x=a:h:b; for i=1:nyp=y(i)+h*feval(fun,x(i),y(i));yc=y(i)+h*feval(fun,x(i+1),yp);y(i+1)=0.5*(yp+yc);end13.% 梯形方法求解一阶常微分方程初值问题% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;% 输出：y--常微分方程在区间[a,b]上各点的数值解；% 调用格式(ex)：y=trapezium(@fun4,0,1,1,10)% Date:2012-1-2function y=trapezium(fun,a,b,y0,n)h=(b-a)/n;y(1)=y0;x=a:h:b;tol=1e-6;for i=1:n%用不动点迭代的方法求解非线性方程：%y(i+1)=y(i)+h*feval(fun,x(i),y(i))/2+h*feval(fun,x(i),y(i+1))/2;iter=1;yy0=1+(i-1)*h;%迭代初始值while iter<500yy1=y(i)+h*feval(fun,x(i),y(i))/2+h*feval(fun,x(i)+h,yy0)/2;if abs(yy1-yy0)<tolbreak;elseyy0=yy1;iter=iter+1;endendif iter>=500printf('not successful!');y=NaN;return;elsey(i+1)=yy1;endend14.% 标准四阶四段龙格-库塔方法求解一阶常微分方程初值问题% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;% 输出：y--常微分方程在区间[a,b]上各点的数值解；% 调用格式(ex)：y=longekuta(@fun4,0,1,1,10)% Date:2012-1-2function y=longekuta(fun,a,b,y0,n)h=(b-a)/n;y(1)=y0;for k=2:n+1x=a+(k-2)*h;k1=h*feval(fun,x,y(k-1));k2=h*feval(fun,x+h/2,y(k-1)+k1/2);k3=h*feval(fun,x+h/2,y(k-1)+k2/2);k4=h*feval(fun,x+h,y(k-1)+k3);y(k)=y(k-1)+(k1+2*k2+2*k3+k4)/6;end15.插值% 拉格朗日插值% 输入：x,y--插值数据点（x，y均为行向量）；xh--要插值的点; % 输出：yh--插值结果；% 调用格式(ex)：yh=lagrange(x,y,xh)% Date:2012-1-2function yh=lagrange(x,y,xh)n=length(x);m=length(xh);yh=zeros(1,m);c1=ones(n-1,1);c2=ones(1,m);for i=1:nxp=x([1:i-1 i+1:n]);yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));%prod对输入的一个向量返回其所有分量的乘积end%拉格朗日调用clear;clc;x=[11 12];y=[2.3979 2.4849];xh=11.75;yh=lagrange(x,y,xh)15_2% 牛顿插值% 输入：x,y--插值数据点（x，y均为行向量）；xh--要插值的点; % 输出：yh--插值结果；% 调用格式(ex)：yh=newtonPol(x,y,xh)% Date:2012-1-2function yh=newtonPol(x,y,xh)n=length(x);p(:,1)=x;p(:,2)=y;for j=3:n+1p(1:n+2-j,j)=diff(p(1:n+3-j,j-1))./(x(j-1:n)-x(1:n+2-j))';%求差商表endq=p(1,2:n+1)';%求牛顿法的系数--取第一行yh=0;m=1;yh=q(1);for i=2:nm=q(i);for j=2:im=m*(xh-x(j-1));%求牛顿法中各多项式值(xh-x0)…(xh-x n-1) endyh=yh+m;%求和end%牛顿插值调用clear;clc;x=[11 12 13];y=[2.3979 2.4849 2.5649];xh=11.75;yh=newtonPol(x,y,xh)。
一些经典的数值分析(matlab程序)
1、牛顿迭代法此方法一般用来求函数的根，速度比较快，程序比较简单。
文件newton.m 内容如下%n表示迭代次数，ret为返回值function ret=newton(n)format longy=inline('x^3+10*x-20');z=inline('3*x^2+10');x0=1.5;for i=1:na=y(x0);b=z(x0);x1=x0-a/b;x0=x1;endret=x12、复合Simpson求积分很简单的一种求定积分的方法，将求积区间分成很多小的曲边梯形，再累加。
文件mulsimpson.m内容如下%复化simpson求积分%a,b表示区间，n表示区间数function ret=mulsimpson(a,b,n)h=(b-a)/n;detsum=0;for i=1:n-1xk=a+i*h;detsum=detsum+fun(xk);endret=h*(fun(a)+fun(b)+2*detsum)/2;%内建函数function z=fun(x)z=exp(-x*x);3、4阶龙格-库塔法求积分此方法速度较快，编程较方便文件step4Runge.m内容如下%4阶Runge-Kutta 法%a,b为积分区间，N为划分数目，y0为初值，函数由fun定义function ret=step4Runge(a,b,N,y0)format longh=(b-a)/N;n=1;x0=a;for n=1:Nx=x0+h;k1=fun(x0,y0);k2=fun(x0+h/2,y0+h*k1/2);k3=fun(x0+h/2,y0+h*k2/2);k4=fun(x0+h,y0+h*k3);y=y0+h*(k1+2*(k2+k3)+k4)/6;x0=xy0=yend%积分函数定义function z=fun(x,y)z=1+y*y;4、Gauss_Seidel迭代解线性方程相比消元法，编程较为容易文件Gauss_Seidel.m内容如下%此函数演示高斯-赛德尔迭代%a表示系数矩阵，b表示值矩阵，n表示系数矩阵阶数，M表示迭代次数%注意b为列向量function y=Gauss_Seidel(a,b,n,M)format longx0=[0;0;0];for k=1:Mfor i=1:ns=0;t=x0(i);for j=1:nif j~=is=s+a(i,j)*x0(j);endendx0(i)=(b(i)-s)/a(i,i);endendy=x0;5、高斯列主消元法此方法为解线性方程组常用方法，先化简增广矩阵，然后回代求解。
追赶法求解三对角线性方程组
追赶法求解三对角线性方程组一 实验目的利用编程方法实现追赶法求解三对角线性方程组。
二 实验内容1、 学习和理解追赶法求解三对角线性方程组的原理及方法；2、 利用MATLAB 编程实现追赶法；3、 举例进行求解，并对结果进行分。
三 实验原理设n 元线性方程组Ax=d 的系数矩阵A 为非奇异的三对角矩阵11222=(1)(n 1)()()a c b a c A a n c b n a n ⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥--⎢⎥⎢⎥⎣⎦………… 这种方程组称为三对角线性方程组。
显然，A 是上下半宽带都是1的带状矩阵。
设A 的前n-1个顺序主子式都不为零，根据定理2.5的推论，A 有唯一的Crout 分解，并且是保留带宽的。
其中L 是下三角矩阵，U 是单位上三角矩阵。
利用矩阵相乘法，可以1112212(1)1u(n 1)()()1l u m l u A LU l n m n l n ⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥==⨯⎢⎥⎢⎥--⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦……………得到：由上列各式可以得到L 和U 。
引入中间量y ，令yUx =，则有：已知L 和d ，可求得y 。
则可得到y 的求解表达式：11/12,3,,()(1)*y()=()[()(1)]/y d l i nm i y i li i di y i di m i y i li==-+=--…1111111/1(2)(1)(1)u (1)(11)/(1)(1)(1)l a l u c u c l mi bi i n a i m i i l i i n ci li ui ui ci li l i a i b i ui=*===≤≤+=+++≤≤-=•=+=+-+Ax LUx Ly d Ly d ====1112222(1)(n 1)(n 1)()()(n)(n)l y d m l y d l n y d m n l n y d ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦……………由y Ux =得：111112221u(n 1)(n 1)(n 1)1(n)(n)u x y u x y x y x y ⎡⎤⎡⎤⎡⎤⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⨯=⎢⎥⎢⎥⎢⎥---⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦………… 可得到X 的求解表达式：()()1,2,,1()()u()(1)x n y n i n n x i y i i x i ==--=-+… 从而得到Ax=d 的解x 。