1、说明：本报告中所有程序均在 Matlab R2006a 中调试，且设置数据显示精度为 16位（format long） 。一、线性方程组求解对于线性方程组 1234123441xx1.用直接法求解；2.用 Jacobi 迭代法求解；3.分别取 ，用 SOR 方法求解比较迭代结果（与精确解比较） 。0.75,1.2,解：1、用直接法求解算法：Gauss 列主元消去法是在 Gauss 消去法中增加选主元的过程，即在第 k 步（k=1,2,3,）消元时，首先在第 k 列主对角元以下（含对角元）元素中挑选绝对值最大的数（即为列主元） ，并通过初等行变换，使得该数位于主对角线上，然后再继续消元。程序：f
2、unction x=gauss(A,b,n)a=A,b;/*a为增广阵%消去过程for k=1:n-1%选主元*/c=0;for q=k:nif abs(a(q,k)cc=a(q,k);l=q;endend%如果主元为0，则矩阵A不可逆if abs(c) A=-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4 b=1 1 1 1n=4gauss(A,b,n)计算结果： 方程组的解为x =-1 -1 -1 -1若采用 LU 分解时，命令如下： A=-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4L U=lu(A) b=1 1 1 1 x=inv(U)*i
3、nv(L)*b计算结果为：x =-1.00000000000000-1.00000000000000-1.00000000000000-1.000000000000002、Jacobi 迭代法算法：对于线性方程组 Ax=b，如果 A 为非奇异方程，则可将 A 分解为：A=D-L-U 其中D 为对角阵，其元素为 A 的对角元素，L 与 U 为 A 的下三角阵和上三角阵。于是Ax=b 化为： ，其中 ， 。11()kkxxDb1()JBL1fDbmatlab 程序：function x=jacobi(A,b,x0)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);
4、B=D(L+U);f=Db;x=B*x0+f;n=1;while norm(x-x0,2)=1.0e-10x0=x;x=B*x0+f;n=n+1;endfprintf(迭代次数为：)nfprintf(n)fprintf(方程组的解为：) A=-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4 b=1 1 1 1 x0=0 0 0 0计算结果：迭代次数为：n =79方程组的解为：ans =-0.99999999986515-0.99999999986515-0.99999999986515-0.999999999865153、SOR 法算法：对于线性方程组Ax=b，如果A为非
5、奇异方程，则可将A 分解为：A=D-L-U其中D 为对角阵，其元素为A的对角元素， L与U 为A 的下三角阵和上三角阵。SOR法是在Gauss-Seidel迭代法的基础上引入松弛因子w ，于是Ax=b 化为： ，1kwkxLf其中 ，1()(wLDwU1()fDwLb程序：function x=sor(A,b,x0,w)D=diag(diag(A);U=-triu(A,1);L=-tril(A,-1);Q=(D-w*L)(1-w)*D+w*U);f=(D-w*L)b*w;x=Q*x0+f;n=1;while norm(x-x0)=1.0e-10x0=x;x=Q*x0+f;n=n+1;endfp
6、rintf(迭代次数为：)nfprintf(n)fprintf(方程组的解为：) A=-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4 b=1 1 1 1 x0=0 0 0 0 sor(A,b,x0,0.75) sor(A,b,x0,1.0) sor(A,b,x0,1.25) sor(A,b,x0,1.5)w=0.75 时计算结果为：迭代次数为：n =74方程组的解为：x =-0.99999999988656-0.99999999989506-0.99999999990292-0.99999999991020w=1.0 时计算结果为：迭代次数为n =42方程组的解为x =
7、-0.99999999992282-0.99999999993294-0.99999999994173-0.99999999994937w=1.25 时计算结果为：迭代次数为n =21方程组的解为x =-1.00000000000666-0.99999999999062-1.00000000000561-0.99999999999973w=1.5 时计算结果为：迭代次数为n =39方程组的解为x =-1.00000000000227-1.00000000001380-0.99999999998230-1.00000000001033方程组的精确解为 x=1 1 1 1,0 x=0:1:10;
8、y=0 0.79 1.53 2.19 2.71 3.03 3.27 2.89 3.06 3.19 3.29;pp=csape(x,y,complet,0.8 0.2)；breaks,coefs,nploys,ncofs,dim=unmkpp(pp) %显示插值节点、系数等计算结果：breaks =0 1 2 3 4 5 6 7 8 9 10coefs =-0.00851403685004 -0.00148596314996 0.80000000000000 0-0.00445788944989 -0.02702807370007 0.77148596314996 0.7900000000000
9、0-0.00365440535039 -0.04040174204975 0.70405614740014 1.53000000000000-0.04092448914854 -0.05136495810093 0.61228944724946 2.190000000000000.10735236194454 -0.17413842554654 0.38678606360200 2.71000000000000-0.26848495862962 0.14791866028708 0.36056629834254 3.030000000000000.42658747257395 -0.65753
10、621560179 -0.14905125697216 3.27000000000000-0.26786493166618 0.62222620212007 -0.18436127045388 2.890000000000000.05487225409078 -0.18136859287848 0.25649633878770 3.060000000000000.05837591530307 -0.01675183060615 0.05837591530307 3.19000000000000nploys =10ncofs =4dim =1注：上述结果的每一行从第一列至第四列分别为 ， ， ,
11、 的系数。3x210xfnplt(pp)得到的三次样条插值曲线如下页图所示：五、解微分方程23,01()1utt1、线性二步法程序：function A=erbufa()h=0.05;u(1)=1;t(1)=0;u(2)=11/9*exp(-0.15)+2/3*(0.05-1/3); %此初值根据精确值计算得到for n=1:19t(n)=h*(n-1);t(n+1)=h*n;t(n+2)=h*(n+1);m=2*h*(2*t(n+1)-3*u(n+1)/3;p=5*h*t(n+2)/6;q=h*(2*t(n)-3*u(n)/6;u(n+2)=(m+p-q+u(n+1)/(1+1.25*h);
12、enda=t u;for k=1:21x(k)=0.05*(k-1);y(k)=11/9*exp(-3*x(k)+2/3*(x(k)-1/3);endB=x y;A=a y;plot(t,u,b,x,y,r);legend(近似值 ,精确值 )title(线性二步法近似值与精确值比较 );计算结果：A =0 1.00000000000000 1.000000000000000.05000000000000 0.86308752674174 0.863087526741740.10000000000000 0.76167414029888 0.749888936388770.1500000000
13、0000 0.67686329201964 0.657101074204390.20000000000000 0.60695104281583 0.581880888559370.25000000000000 0.55004943137387 0.521781342239020.30000000000000 0.50451915072020 0.474696250794070.35000000000000 0.46892718879924 0.438812804469190.40000000000000 0.44202112817630 0.412570703448250.4500000000
14、0000 0.42270637968187 0.394626985233870.50000000000000 0.41002629325624 0.383825751292530.55000000000000 0.40314477498604 0.379172110536480.60000000000000 0.40133109473146 0.379809752270830.65000000000000 0.40394660828200 0.385001643050180.70000000000000 0.41043315277373 0.394113412309200.7500000000
15、0000 0.42030290450517 0.406599052242280.80000000000000 0.43312951486180 0.421988609575950.85000000000000 0.44854036328306 0.439877591779190.90000000000000 0.46620978650318 0.459917848904140.95000000000000 0.48585316103681 0.481809725513691.00000000000000 0.50722173138419 0.50529530578294注：上述结果中第一列为
16、t 的值，第二列为 采用线性二步法法计算出的近似解，第三列为精确解。2、经典 Runge-Kutta 方法程序：function A=rungekutta() h=0.05;u(1)=0;t(1)=0;for n=1:20t(n)=h*(n-1);k1=2*t(n)-3*u(n);k2=2*(t(n)+0.5*h)-3*(u(n)+0.5*h*k1);k3=2*(t(n)+0.5*h)-3*(u(n)+0.5*h*k2);k4=2*(t(n)+h)-3*(u(n)+h*k3);u(n+1)=u(n)+h/6*(k1+2*k2+2*k3+k4);enda=t u;for k=1:21x(k)=0
17、.05*(k-1);y(k)=11/9*exp(-3*x(k)+2/3*(x(k)-1/3);endB=x y;A=a y;plot(t,u,b,x,y,r);legend(近似值 ,精确值)title(Runge-Kutta方法近似值与精确值比较); A=rungekutta()计算结果：A=0 1.00000000000000 1.000000000000000.05000000000000 0.86308828125000 0.863087526741740.10000000000000 0.74989023521179 0.749888936388770.15000000000000
18、0.65710275106600 0.657101074204390.20000000000000 0.58188281294427 0.581880888559370.25000000000000 0.52178341265656 0.521781342239020.30000000000000 0.47469838922470 0.474696250794070.35000000000000 0.43881495179498 0.438812804469190.40000000000000 0.41257281570093 0.412570703448250.45000000000000
19、0.39462903052143 0.394626985233870.50000000000000 0.38382770728803 0.383825751292530.55000000000000 0.37917396243216 0.379172110536480.60000000000000 0.37981149111660 0.379809752270830.65000000000000 0.38500326440906 0.385001643050180.70000000000000 0.39411491517368 0.394113412309200.75000000000000
20、0.40660043816504 0.406599052242280.80000000000000 0.42198988197617 0.421988609575950.85000000000000 0.43987875539243 0.439877591779190.90000000000000 0.45991890934932 0.459917848904140.95000000000000 0.48181068895509 0.481809725513691.00000000000000 0.50529617866925 0.50529530578294注：上述结果中第一列为 t 的值，第二列为 Runge-Kutta方法计算出的近似解，第三列为精确解。从上述计算结果中可以看出四阶计算得到的近似解与精确解的误差很小，小于 10-5，在图形上已经无法显示出误差。精确解：3121()9tuet01tk=1:21x(k)=0.05*(k-1);y(k)=11/9*exp(-3*x(k)+2/3*(x(k)-1/3);plot(x,y,r);title(精确解)
