数值分析——编程作业.docx

数值分析实验报告第二章： 解线性方程组的直接方法2、试用MATLAB软件编程实现追赶法求解三对角方程组的算法，并考虑梯形电阻电路问题，电路如下：其中电路中的各个电流，须满足下列线性方程组：设，运用求各段电路的电流量。解：上述方程组可用矩阵表示为：MatLab程序： %赋初值； 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; 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:1 d(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即 。1、试分别用（1）Jacobi迭代法；（2）Gauss-Seidel迭代法；（3）共轭梯度法解线性方程组迭代初始向量取=0，0，0，0，0。解：实验步骤及程序、结果 取要求达到的精度。以下程序中的均表示迭代次数。（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:5 sum=0;for j=1:5 if 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 & k1e-8 & k=deg & k=deg & kN k=k+1; eigenvalue0=eigenvalue;xr,r=max(abs(x); a=xr; y=x/a; x=B*y; eigenvalue=xr; end; eigenvalue x k输出结果：eigenvalue = 0.002637740424370x = -0.001704294528000 0.000852172462123 0.008922284278451 0.000568148573659k = 68第五章 插值法3、已知直升机旋转机翼外形轮廓线上某些型值点的数据： 1 520 02 280 -303 156.6 -364 78 -355 39.62 -28.446 3.1 -9.47 0 08 3.1 9.49 39.62 28.4410 78 3511 156.6 3612 280 3013 520 0由表中13个节点，用样条插值的方法，增加平面点为101各点，绘制出较平滑的机翼外形曲线图。解：了解三次样条插值方法的使用。程序如下：x=-520,-280,-156.6,-78,-39.62,-3.1,0,3.1,39.62,78,156.6,280,520;y=0,-30,-36,-35,-28.44,-9.4,0,9.4,28.44,35,36,30,0;n=13;%求解Mfor i=1:1:n-1h(i)=x(i+1)-x(i);endfor i=2:1:n-1a(i)=h(i-1)/(h(i-1)+h(i);b(i)=1-a(i);c(i)=6*(y(i+1)-y(i)/h(i)-(y(i)-y(i-1)/h(i-1)/(h(i-1)+h(i);enda(n)=h(n-1)/(h(1)+h(n-1);b(n)=h(1)/(h(1)+h(n-1);c(n)=6/(h(1)+h(n-1)*(y(2)-y(1)/h(1)-(y(n)-y(n-1)/h(n-1);A(1,1)=2;A(1,2)=b(2);A(1,n-1)=a(2);A(n-1,n-2)=a(n);A(n-1,n-1)=2;A(n-1,1)=b(n);for i=2:1:n-2A(i,i)=2;A(i,i+1)=b(i+1);A(i,i-1)=a(i+1);endC=c(2:n);C=C;m=AC;M(1)=m(n-1);M(2:n)=m;xx=-520:10.4:520;for i=1:51 for j=1:1:n-1 if x(j)=xx(i) & xx(i)3*dx Fx=0; T0=T; for k=1:2(m-1) Fx=Fx+xf(a+(2*k-1)*h); end T=T0/2+h*Fx; m=m+1; h=h/2; end x(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=Fy+yf(a+(2*k-1)*h); end T=T0/2+h*Fy; n=n+1; h=h/2; end y(j)=T; j=j+1;endplot(x,y,o-g);在命令窗口中输入x,y=fuhuatixing()可得x与y的值，即曲线上点的坐标，在本例中s从-5到5步长为0.1计算了101个点的坐标值，得到曲线图形如下图所示：第八章 非线性方程及非线性方程组的解法1. 求下列方程的非零根分析：本题拟采用牛顿法对该非线性方程进行求解，而本函数f(x)为奇函数，为了使得牛顿迭代公式初始值比较靠近该方程的非零根，应该确定该方程非零根的大概范围。故首先用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函数图如下：由上图可以很明显看出，原方程有一个零根，两个关于原点对称的非零实根。由图估测原方程的正实根的x坐标值在770附近。并且经验证f(770)= 1.07710，而f(760)= -1.00550。所以正实根落在区间760,770上。现取初始值=765。误差限errorlim=，最大迭代次数N=100M文件代码如下：%子函数非线性函数f3 function y=f3(x)y=log(513+0.6651*x)/(513-0.6651*x)-x/(1400*0.0918);end%子函数非线性函数一阶导数df3 function y=df3()syms x1y=log(513+0.6651*x1)/(513-0.6651*x1)-x1/(1400*0.0918);y=diff(y);end%主程序-牛顿法求解非线性方程clearclcformat longx0=765;N=100;errorlim=10(-5);x=x0-f3(x0)/subs(df3(),x0);n=1; while nerrorlim n=n+1; else break; end x0=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、 设常微分方程初值问题其精确解为。选取步长h使四阶Adams预测-校正算法和经典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=n n=round;endfor i=1:n k1=fxy(x(i),yrk4(i); k2=fxy(x(i)+h/2,yrk4(i)+h*k1/2); k3=fxy(x(i)+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 on%绘图-数值解plot(x1,yrk4,xr)grid onlegend(精确解,数值解)title(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:n p(i+1)=yadms(i)+(h/24)*(55*fxy(x(i),yadms(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)f=-y+2*cos(x);%主程序-Adams预测校正算法clearclca=0;b=pi;N=50;y0=1;yrk4