数值分析计算方法实验报告.doc

课程名称计算方法实验项目名 称线性方程组的数值解法实验项目类型验证演示综合设计其他指导教师胡小兵成 绩第1页，共14页1实验目的：（1）高斯列主元消去法求解线性方程组的过程（2）熟悉用迭代法求解线性方程组的过程（3）设计出相应的算法，编制相应的函数子程序2实验内容分别用高斯列主元消去法 ,Jacobi迭代法,Gauss-Saidel迭代法,超松弛迭代法求解线性方程组3实验环境及实验文件存档名写出实验环境及实验文件存档名4.实验结果及分析输出计算结果，结果分析和小结等。解：1. 高斯列主元消去法： %用高斯列主元消去法解实验1%高斯列主元消元法求解线性方程组Ax=b%A为输入矩阵系数,b为方程组右端系数%方程组的解保存在x变量中format long;A=2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13b=10,5,-2,7m,n=size(A);if m=n error(矩阵A的行数和列数必须相同); return;endif m=size(b) error(b的大小必须和A的行数或A的列数相同); return;endif rank(A)=rank(A,b) error(A矩阵的秩和增广矩阵的秩不相同,方程不存在唯一解); return;endc=n+1;A(:,c)=b; for k=1:n-1r,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,:); end A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c); endendx=zeros(length(b),1);x(n)=A(n,c)/A(n,n);for k=n-1:-1:1 x(k)=(A(k,c)-A(k,k+1:n)*x(k+1:n)/A(k,k);enddisp(X=);disp(x);format short;输出结果： 2Jacobi迭代法：%Jacobi迭代法求解实验1% A为方程组的增广矩阵clc;A=2 10 0 -3 10;-3 -4 -12 13 5;1 2 3 -4 -2;4 14 9 -13 7MAXTIME=50;eps=1e-5;n,m=size(A); x=zeros(n,1); y=zeros(n,1); k=0; disp(迭代过程X的值情况如下:) disp(X=); while 1 disp(x); for i=1:1:n s=0.0; for j=1:1:n if j=i s=s+A(i,j)*x(j); end y(i)=(A(i,n+1)-s)/A(i,i); end end for i=1:1:n maxeps=max(0,abs(x(i)-y(i); end if maxepsMAXTIME error(超过最大迭代次数,退出); return; end end输出结果：由于不收敛，故出现上述情况。3.Gauss-Saidel迭代法%Gauss_Seidel迭代法求解实验1% A为方程组的增广矩阵clc;format long;A=2 10 0 -3 10;-3 -4 -12 13 5;1 2 3 -4 -2;4 14 9 -13 7n,m=size(A);Maxtime=50;Eps=10E-5;x=zeros(1,n);disp(x=);for k=1:Maxtime disp(x); for i=1:n s=0.0; for j=1:n if i=j s=s+A(i,j)*x(j); end end x(i)=(A(i,n+1)-s)/A(i,i); end if sum(x-floor(x).2)Eps break; end;end;X=x;disp(迭代结果:);Xformat short;输出结果：因为不收敛，故出现上述情况。4.超松弛迭代法：%SOR法求解实验1%w=1.45%方程组系数矩阵clc;A=2,10,0,-3;-3,-4,-12,13;1,2,3,-4;4,14,9,-13b=10,5,-2,7w=1.45;Maxtime=100;Eps=1E-5;format long;n=length(A);k=0;x=ones(n,1);y=x;disp(迭代过程:);disp(x=);while 1 y=x; disp(x); for i=1:n s=b(i); for j=1:n if j=i s=s-A(i,j)*x(j); end end if abs(A(i,i)=Maxtime error(已达最大迭代次数或矩阵系数近似为0,无法进行迭代); return; end s=s/A(i,i); x(i)=(1-w)*x(i)+w*s; end if norm(y-x,inf)Eps break; end k=k+1;enddisp(最后迭代结果:);X=xformat short；输出结果：由于不收敛，故出现上述情况。课程名称计算方法实验项目名 称插值方法实验项目类型验证演示综合设计其他指导教师胡小兵成 绩第7页，共14页1.实验目的：（1） 学会拉格朗日插值、牛顿插值等基本方法（2） 设计出相应的算法，编制相应的函数子程序（3） 会用这些函数解决实际问题2实验内容（1）设计拉格朗日插值算法，编制并调试相应的函数子程序（2）设计牛顿插值算法，编制并调试相应的函数子程序（3）给定函数四个点的数据如下：X1.12.33.95.1Y3.8874.2764.6512.117 试用拉格朗日插值确定函数在x=2.101，4.234处的函数值。4）已知用牛顿插值公式求的近似值。3实验原理写出本次实验所用算法的算法步骤叙述或画出算法程序框图4实验环境及实验文件存档名写出实验环境及实验文件存档名5.实验结果及分析输出计算结果，结果分析和小结等。（1）拉格朗日法程序：function lagrange(x)format long;x0=1.1 2.3 3.9 5.1;y0=3.887 4.276 4.651 2.117;n=length(x0);s=0;for j=0:(n-1) t=1; for i=0:(n-1) if i=j t=t*(x-x0(i+1)/(x0(j+1)-x0(i+1); end end s=s+t*y0(j+1);ends format short;运行结果：（2）牛顿差值法：function Newton(x)format long;%显示15位x0=1 4 9;%x的值y0=1 2 3;%y的值x=5;%插值点n=max(size(x0);y=y0(1);%迭代初始值disp(y);s=1;dx=y0;for i=1:n-1%构造差商表 dx0=dx; for j=1:n-i dx(j)=(dx0(j+1)-dx0(j)/(x0(i+j)-x0(j); end df=dx(1); s=s*(x-x0(i); y=y+s*df;%计算 disp(y);end运行结果：实际结果是2.236066797749979.有一定的误差。总体来说还是很准确的课程名称计算方法实验项目名 称数值微积分实验项目类型验证演示综合设计其他指导教师胡小兵成 绩第3页，共14页1实验目的：（1）学会复化梯形、复化辛浦生求积公式的应用（3）设计出相应的算法，编制相应的函数子程序（4）会用这些函数解决实际问题2实验内容（1）设计复化梯形公式求积算法，编制并调试相应的函数子程序（2）设计复化辛浦生求积算法，编制并调试相应的函数子程序（4）分别用复化梯形公式和复化辛浦生公式计算定积分取n=2，4，8，16，精确解为0.94608313、 实验原理写出本次实验所用算法的算法步骤叙述或画出算法程序框图4实验环境及实验文件存档名写出实验环境及实验文件存档名5.实验结果及分析输出计算结果，结果分析和小结等。本实验是用MATLAB来完成。梯形法存档名为trap.m 辛浦生法存档名为simpson.m复化梯形公式求积算法：function T=trap(n)f=sym(sin(x)/x);a=0;b=1;h=(b-a)/n;T=0;for k=1:(n-1) x0=a+h*k; T=T+limit(f,x0);endT=h*(limit(f,a)+limit(f,b)/2+h*T;T=double(T);复化辛浦生法：function S=simpson(n)f=sym(sin(x)/x);a=0;b=1;h=(b-a)/(2*n);s1=0;s2=0;for k=1:n x0=a+h*(2*k-1); s1=s1+limit(f,x0);endfor k=1:(n-1) x0=a+h*2*k; s2=s2+limit(f,x0);endS=h*(limit(f,a)+limit(f,b)+4*s1+2*s2)/3;S=double(S);运行结果 课程名称计算方法实验项目名 称常微分方程的数值解法实验项目类型验证演示综合设计其他指导教师胡小兵成 绩第12页1实验目的：（1）学会四阶龙格-库塔方法的使用（2）设计出相应的算法，编制相应的函数子程序（3）会用这些函数解决实际问题2实验内容（1）分别取h=0.05，N=10；h=0.025，N=20；h=0.01，N=50，用四阶龙格-库塔方法求解微分方程初值 问题：y=-50y，y(0)=10（2）某跳伞者在t=0时刻从飞机上跳出，假设初始时刻的垂直速度为0，且跳伞者垂直下落。已知空气阻力为F=cv2，其中c为常数，v为垂直速度，向下方方向为正。写出此跳伞者的速度满足的微分方程；若此跳伞者的质量为M=70kg，且已知c=0.27kg/m，利用四阶龙格-库塔公式计算t=20s的速度（取h=0.1s）3.实验原理写出本次实验所用算法的算法步骤叙述或画出算法程序框图3、 实验结果及分析输出计算结果，结果分析和小结等。（1）程序：function RK(h,n)F=-50*y;a=0;b=a+h*n;X1=a:h:b;Y1=zeros(1,n+1);Y1(1)=10;for i=1:n x=X1(i); y=Y1(i); K1=h*eval(F); x=x+h/2; y=y+K1/2; K2=h*eval(F); x=x; y=Y1(i)+K2/2; K3=h*eval(F); x=X1(i)+h; y=Y1(i)+K3; K4=h*eval(F); Y1(i+1)=Y1(i)+(K1+2*K2+2*K3+K4)/6;endX1disp(Y1);endh=0.05，N=10时，运行结果：h=0.025,n=20时，运行结果：h=0.01，N=50时，输出结果：（2）微分方程为：v=9.8-0.003857v2,v(0)=0程序：function velosity(h,n)F=9.8-0.003857