1、数值分析实验报告第二章： 解线性方程组的直接方法2、试用 MATLAB 软件编程实现追赶法求解三对角方程组的算法，并考虑梯形电阻电路问题，电路如下：其中电路中的各个电流 ， ， 须满足下列线性方程组：1i28iRVi 2 1 053 42ii2 53 0 64ii 75 2 86ii05 7设 ， ，运用求各段电路的电流量。V2027R解： 148.上述方程组可用矩阵表示为： 00148.5200520052 87654321iiiiMatLab 程序：%赋初值；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;
2、d=8.1481 0 0 0 0 0 0 0;%三对角方程的追赶法for i=2:8 %“追”的过程；a(i)=a(i)/b(i-1);b(i)=b(i)-c(i-1)*a(i);d(i)=d(i)-a(i)*d(i-1);end;d(8)=d(8)/b(8); %“赶”的过程；for i=7:-1:1d(i)=(d(i)-c(i)*d(i+1)/b(i);end;x=d;x程序运行结果：x =8.1477 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477即 。(A) 047.19256.307.41887654321iiiiI1、试分别用（1）
3、Jacobi 迭代法；（2）Gauss-Seidel 迭代法；（3）共轭梯度法解线性方程组12741215342712394042xx迭代初始向量取 = 0，0，0，0，0 。)(xT)解：实验步骤及程序、结果取要求达到的精度 。以下程序中的 均表示迭代次数。81k（1）Jacobi 迭代法 MatLab 源程序。format longA=10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15;b=12,-27,14,-17,12;x0=0,0,0,0,0;x1=x0;Nmax=1000;k=0;for i=1:5sum=0;f
4、or j=1:5if j=isum=sum+A(i,j)*x0(j);end;end;x1(i)=(b(i)-sum)/A(i,i);end;while abs(norm(x1-x0,inf)1e-8 T0=T;for k=1:2(m-1)Fx=Fx+xf(a+(2*k-1)*h);endT=T0/2+h*Fx; m=m+1;h=h/2; endx(i)=T;i=i+1;endj=1;for s=-5:0.1:5 n=1;b=s;a=0;h=(b-a)/2;T=h*(yf(a)+yf(b);T0=5;while abs(T-T0)3*dy Fy=0;T0=T;for k=1:2(n-1)Fy=
5、Fy+yf(a+(2*k-1)*h);endT=T0/2+h*Fy; n=n+1;h=h/2; endy(j)=T;j=j+1;endplot(x,y,o-g);在命令窗口中输入x,y=fuhuatixing()可得 x 与 y 的值，即曲线上点的坐标，在本例中 s 从-5 到 5 步长为 0.1 计算了101 个点的坐标值，得到曲线图形如下图所示：第八章 非线性方程及非线性方程组的解法1. 求下列方程的非零根 5130.6()ln0140.98xfx分析：本题拟采用牛顿法对该非线性方程进行求解，而本函数 f(x)为奇函数，为了使得牛顿迭代公式初始值比较靠近该方程的非零根，应该确定该方程非零根
6、的大概范围。故首先用 Matlab 画图指令画出该函数的曲线图。因注意到 513-0.6651x0，则可知|x| x=-770:770; y=log(513+0.6651*x)./(513-0.6651*x)-x/(1400*0.0918); plot(x,y) grid函数图如下：-800 -600 -400 -200 0 200 400 600 800-3-2-10123由上图可以很明显看出，原方程有一个零根，两个关于原点对称的非零实根。由图估测原方程的正实根的 x 坐标值在 770 附近。并且经验证 f(770)= 1.07710，而 f(760)= -1.0055errorlimn=n
7、+1;elsebreak;endx0=x;enddisp(非线性方程求解完毕！)disp(迭代次数 : n=,num2str(n)disp(所求非零根 : 正根 x1=,num2str(x), 负根 x2=,num2str(-x)运行结果：非线性方程求解完毕！迭代次数: n=5所求非零根: 正根 x1=767.3861 负根 x2=-767.3861结果验证： f3(x)ans =-4.440892098500626e-015 f3(-x)ans =4.440892098500626e-015第九章 常微分方程数值解法1、 设常微分方程初值问题 2cos,0()1yx其精确解为 。选取步长 h
8、 使四阶 Adams 预测-校正算法和经cosinyx典 RK 法均稳定，分别用这两种方法求解微分方程，将数值解和精确解进行比较，输出结果。其中多步法需要的初值由经典 RK 法提供。解：M 文件代码：%子函数-四阶 RK 算法function yrk4=rk4(a,b,y0,n)h=(b-a)/n;x=a:h:b;yrk4=zeros(n+1,1);yrk4(1)=y0;%round控制 RK-4 算法计算步数if round=nn=round;endfor i=1:nk1=fxy(x(i),yrk4(i);k2=fxy(x(i)+h/2,yrk4(i)+h*k1/2);k3=fxy(x(i)
9、+h/2,yrk4(i)+h*k2/2);k4=fxy(x(i)+h,yrk4(i)+h*k3);yrk4(i+1)=yrk4(i)+h*(k1+2*k2+2*k3+k4)/6;end%子函数 -fxy 函数function f=fxy(x,y)f=-y+2*cos(x);%主程序 -经典 RK 法求解常微分方程clearclca=0;b=pi;N=50;n=1;y0=1;yrk4=rk4(a,b,y0,N,N);%绘图 -精确解y=dsolve(Dy=-y+2*cos(x),y(0)=1,x);x1=0:(pi/N):pi;y1=subs(y,x1);plot(x1,y1,-b)hold o
10、n%绘图 -数值解plot(x1,yrk4,xr)grid onlegend(精确解 ,数值解)title(RK-4 数值解与精确解)运行结果：0 0.5 1 1.5 2 2.5 3 3.5-1-0.500.511.5 RK-4与与与与与与与与与与与与与%子函数 -Adams 预测校正算法function yadms=adms(a,b,y,n)h=(b-a)/n;x=a:h:b;yadms=y;zeros(n-3,1);N=length(yadms);p=zeros(N,1);m=p;c=m;for i=4:np(i+1)=yadms(i)+(h/24)*(55*fxy(x(i),yadms(
11、i)-59*fxy(x(i-1),.yadms(i-1)+37*fxy(x(i-2),yadms(i-2)-9*fxy(x(i-3),yadms(i-3);m(i+1)=p(i+1)+(251/270)*(c(i)-p(i);c(i+1)=yadms(i)+(h/24)*(9*fxy(x(i+1),m(i+1)+19*fxy(x(i),.yadms(i)-5*fxy(x(i-1),yadms(i-1)+fxy(x(i-2),yadms(i-2);yadms(i+1)=c(i+1)-(19/270)*(c(i+1)-p(i+1);end%子函数 -fxy 函数function f=fxy(x,y
12、)f=-y+2*cos(x);%主程序 -Adams 预测校正算法clearclca=0;b=pi;N=50;y0=1;yrk4=rk4(a,b,y0,N,4);yadms=adms(a,b,yrk4(1:4),N);%绘图 -精确解y=dsolve(Dy=-y+2*cos(x),y(0)=1,x);x1=0:(pi/N):pi;y1=subs(y,x1);plot(x1,y1,-b)hold on%绘图 -数值解plot(x1,yadms,xr)grid onlegend(精确解 ,数值解)title(Adams 预测校正数值解与精确解 )运行结果：0 0.5 1 1.5 2 2.5 3 3.5-1.5-1-0.500.511.5 Adams与与与与与与与与与与与与与与与与与由实验结果可以看出，该 2 种算法是稳定的。求解效果良好。
