数值分析实验报告(同名6593).docx

《数值分析实验报告(同名6593).docx》由会员分享，可在线阅读，更多相关《数值分析实验报告(同名6593).docx（28页珍藏版）》请在三一办公上搜索。

1、学生实验报告实验课程名称开课实验室数学与统计学院实验室学 院 2010年级数学与应用数学专业班01班学生姓名 学号开课时间 2012至 2013学年第一学期总成绩教师签名程数值分析实验项目Gauss消元法实验项目类型名称名 称验证演示综合设计其他指导教师何光辉成绩V一、实验目的：（1）高斯列主元消去法求解线性方程组的过程（2）熟悉用迭代法求解线性方程组的过程（3）设计出相应的算法，编制相应的函数子程序二、实验内容分别用高斯列主元消元法和直接消元法求解线性方程组:_ 2100-3_x1 10 -3-4-1213x 2=5123-4x3-24149-13x1- 47三、实验原理对于线性方程组a x

2、 + a x + + a x =b11 112 21nn1a x + a x + + a x =b 21 122 22 nn2( 1)a x + a x + + a x =bn1 1 n 2 2nn n n常记为矩阵形式Ax=b（2）根据高等代数的知识，若| A|。0，上式的解存在且唯一。（1） Gauss直接消元法考虑上述线性方程组的增广矩阵A: b，对增广矩阵进行行变换，将（2）式化为等价的三角形方阵，然后回代解之，这就是Gauss消元法。具体如下: a）消元 令a（i）= a , i, j = 1,2，n ； b（i）= b, i = 1,2，n 对k= 1到n-1，若a。0，进行a (

3、k)i = k + 1,k + 2,，nl =-ik, a (k) kk a (k+1) ika (k+1)=0,=a(k) 一l xa(k),i = k +1,k + 2,，ni,j = k + 1,k + 2,，njb( k+1)iijikkj=b(k) 一 l x bk ,iikki = k + 1,k + 2,，nb)回代，若a如。0_ bnx 一n bn/nnx = (b(i)-不 a( i)x )i ai iij jiij=k+l(2)Gauss列主元消元法设列主元消元法已完成Ax = b的第k-1 (1 k n -1)次消元，的到方程组Ax = b A（ k）x = b（ k）在

4、进行第k次消元前，先进行2个步骤:a）在ak)至a(?这一列内选出最大值，即a(k)=maxik,kk i na（ k）ik若 a(k)k = 0，此时 A = 0方程组无确定解，应给出退出信息。b）若a(k)。0，则交换第i行和k行，然后用Gauss消元法进行消元。k四、MATLAB软件实现(1) 写出Gauss消元法和列主元消元法实现的MATLAB函数根据以上的算法，写出如下程序：%Gauss 消元法 %function y=Gauss1(A,b)m,n=size(A);%检查系数正确性if m=nerror，矩阵A的行数和列数必须相同)；return;endif m=size(b)err

5、or b的大小必须和A的行数或A的列数相同);return;end%再检查方程是否存在唯一解if rank(A)=rank(A,b)error A矩阵的秩和增广矩阵的秩不相同，方程不存在唯一解)； return;end%这里采用增广矩阵行变换的方式求解c=n+1;A(:,c)=b;%消元过程for k=1:n-1A(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c); end%回代结果-WORD格式-可编辑-专业资料- x=zeros(length(b),1);x(n)=A(n,c)/A(n,n); for k=n-1：-1:1x(k

6、) = (A(k,c)-A(k,k+1:n)*x(k+1:n)/A(k,k); end%显示计算结果 %disp(x=); %disp(x);y=x;% %高斯列主元消元法求解线性方程组Ax=b% %A为输入矩阵系数，b为方程组右端系数 %方程组的解保存在x变量中 function y=Gauss_line(A,b)format long; %设置为长格式显示，显示15位小数 m,n=size(A);%先检查系数正确性 if m=nerror，矩阵A的行数和列数必须相同)； return;end if m=size(b)error b的大小必须和A的行数或A的列数相同)； return;end

7、%再检查方程是否存在唯一解 if rank(A)=rank(A,b)error A矩阵的秩和增广矩阵的秩不相同，方程不存在唯一解)； return;endc=n+1;A(:,c)=b; %(增广) for k=1:n-1 r,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,:) %换行 endA(k+1:n, k:c)=A(k+1:n, k:c)-(A(k+1:n,k)/ A(k,k)*A(k, k:c); %消去 end end x=zeros(length(b),1); % 回代求

8、解军 x(n)=A(n,c)/A(n,n); for k=n-1:-1:1x(k) = (A(k,c)-A(k,k+1:n)*x(k+1:n)/A(k,k);endy=x;format short; %设置为默认格式显示，显示5位(2) 建立MATLAB界面利用MATLAB的GUI建立如下界面求解线性方程组：详见程序。五、计算实例、数据、结果、分析下面我们对以上的结果进行测试，求解:一 2100-3_X1 10 -3-4-1213X 2=5123-4X3-24149-13X1- 47输入数据后点击I心好履I和G-L好方豆，得到如下结果：Gauss直接法求解结果122 L*1Gauss列主元法求

9、解结果Copyright: yfg Td: 13500398935更改以上数据进行测试，求解如下方程组：得到如下结果:4 3 2 13 4 3 22 3 4 312 3 4尤1一1X12=X一13X-11- 4 ,求解线性方程组Ax=b清按matlab中输入矩阵方式输入矩阵A4321;3432;2343;1 234请按matlab中输入矩阵方式输入bGauss直接法求蟀结果G和流刁Gauss列主元法求解结果-G-球犀方程I 1 - -127756e-0171-1-1 33230162.7756&fl1711-1 3323&XJ16Cmvrinht vfn TeIT 苗fin舶舶35六、实验中遇

10、到的问题及解决办法在本实验中，遇到的问题主要有两个：(1) 如何将上述的Gauss消元法的算法在MATLAB中实现针对此问题我借鉴了网上以及 课本上的算法的MATLAB实现的程序；(2) 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完 成了此项任务。七、实验结论通过以上的测试，我们发现以上算法和程序能够求出线性方程组的比较精确解。八、参考文献1杨大地，王开荣.2006.数值分析.北京:科学出版社何光辉.2008.数值分析实验.重庆大学数理学院数学实验教学中心(3) 百度文库，百度知道教师签名年 月

11、日课程名称数值分析实验项目名 称插值方法实验项目类型验证演示综合设计其他指导教师何光辉成绩V、实验目的:（1）（2）（3）实验内容（1）设计拉格朗日插值算法，编制并调试相应的函数子程序（2）设计牛顿插值算法，编制并调试相应的函数子程序（3）给定函数四个点的数据如下：X1.12.33.95.1Y3.8874.2764.6512.117试用拉格朗日插值确定函数在x=2.101, 4.234处的函数值。学会拉格朗日插值、牛顿插值等基本方法 设计出相应的算法，编制相应的函数子程序 会用这些函数解决实际问题（4）已知t1 = 1,4 = 2,9 =3,用牛顿插值公式求5的近似值。三、实验原理(1) 拉格

12、朗日插值n 次拉格朗日插值多项式为：Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+.+ynln(x)n=1 时，称为线性插值，L1(x)=y0(x-x1)/(x0-x1)+ y1(x-x0)/(x1-x0)=y0+(y1-x0)(x-x0)/(x1-x0) n=2时，称为二次插值或抛物线插值，精度相对高些 L2(x)=y0(x-x1)(x-x2)/(x0-x1)/(x0-x2)+y1(x-x0)(x-x2)/(x1-x0)/(x1-x2)+y2(x-x0)(x-x1)/(x2-x0)/(x2-x1) 对节点xi(i=0,1,.,n冲任一点xk(0=k=n)作一 n次多项式lk(x

13、k)，使它在该点上取值为1，而在其 余点 xi(i=0,1,.,k-1,k+1,.,n)上为 0，则插值多项式为Ln(x)=y0l0(x)+y1l1(x)+y2l2(x)+.+ynln(x) 上式表明：n 个点 xi(i=0,1,.,k-1,k+1,.,n)都是 lk(x)的零点。(2) 牛顿插值插商公式kf (x )f x , x , x ,.x = 2 j;1 2 3 k (x x ).(x x)(x x).(x x )j=1j 1 jj1jj1j kNewton插值多项式为N (x) = f (x ) + f x , x (x - x ) + f x , x , x (x - x )(x

14、 - x ) +n112112312f x , x , x ,.x (x x )(x x ).(x x);123 n12n1四、MATLAB软件实现（1）分别写出lagrange插值法和Newton插值法的求解函数%lagrange 插值法求解函数 %x, y为初始数据，z为插值点 function z=lagrange(x,y,a) format long; % 显示 15位 n=length(x);% 取长度%初始计算s = 0;%进入公式计算for j=0:(n-1)t=1;for i=0:(n-1) if i=jt=t*(a-x(i+1)/(x(j+1)-x(i+1);endends=

15、s+t*y(j+1);endz=s; %显示输出结果 format short;%Newton 插值法求解函数 %x，y为初始数据，z为插值点 function j=Newton(x,y,z) n=max(size(x);l=1;a=y(1);B=a;s = 1; %一次因子的乘积，预设为1 dx=y; %差商for i=1:n-1dx0=dx;for j=1:n-idx(j) = (dx0(j+1)-dx0(j)/(x(i+j)-x(j);enddf=dx(1);s=s*(z-x(i); %一次因子乘积a=a+s*df; %计算各次Newton插值的值 l=l+1;B=a%结果保存在变量B中

16、endj=B;(2)建立界面利用MATLAB中的GUI编程建立如下界面：Lagrange插值结果4.3017执咛指由1_Newton插值结果ICrnvrinlTt- vtn TeIT 和忡沼节同理可以测试（4）中的忘的值。六、实验中遇到的问题及解决办法在本实验中，遇到的问题主要有两个：（3）如何将上述的插值的算法在MATLAB中实现针对此问题我借鉴了网上以及 课本上的算法的MATLAB实现的程序；（4）如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完 成了此项任务。七、实验结论通过以上的测试，我们发现

17、以上算法和程序能够求出插值的比较精确解。八、参考文献1杨大地，王开荣.2006.数值分析.北京:科学出版社何光辉.2008.数值分析实验.数理学院数学实验教学中心3百度文库，百度知道教师签名年 月 日课程名称数值分析实验项目名 称数值微积分实验项目类型验证演示综合设计其他指导教师何光辉成绩V一、实验目的：（1）学会复化梯形、复化辛浦生求积公式的应用（2）设计出相应的算法，编制相应的函数子程序（3）会用这些函数解决实际问题实验内容（1）设计复化梯形公式求积算法，编制并调试相应的函数子程序（2）设计复化辛浦生求积算法，编制并调试相应的函数子程序（4）分别用复化梯形公式和复化辛浦生公式计算定积分!

18、sin 工,j1 dx三、实验原理（1）复化梯形求积公式h u; x u a; T u 0n图1复化梯形求积公式算法的流程图Step1给出被积函数f （x）、区间a,b端点a,b和等分数n ；h aStep2 求出 Xk = kh, h =；Step3 计算 f (a), f (b),区 f (xk)；k=0f (X ) + f (b) kStep4 得 T =1 h f (a) + 云 n 2k=0(2)复化辛普森求积公式图2复化辛普森求积公式算法的流程图Stepl给出被积函数f (x)、区间a,b端点a,b和等分数n ；.h aStep2 求出 Xk = kh, h =；Step3计算寸（

19、饥，云/（号，习/（:+； k=0k=0k 2Step4 得 S = f (a) + 4云 f (x) + 21 f (x ) + f (b)n 6k+】kk=02k=1四、MATLAB软件实现(1) 分别写出复化梯形和复化辛浦生求积的求解函数 % 复化梯形公式求积分值 % function T=trap(f,a,b)%f为积分函数%a,b为积分区间%n是等分区间份数n=200;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);%Sim

20、pson 公式求积分值 % function S=simpson(f,a,b)%f为积分函数%a,b为积分区间%n是等分区间份数n=200;h=(b-a)/(2*n); % 步长s1=0;s2 = 0;for k=1:nx0=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);(2)建立界面利用MATLAB中的GUI编程建立如下界面：积分求解请输入非常数的积分函数请输入积分区间一Simps

21、on法求解结果Simpsicnt复化梯形公式求解结果|便化盘移详见程序。五、计算实例、数据、结果、分析下面我们对以上的问题进行测试：输入数据：请输入非常数的积分函数请输入积分区间|。计算结果如下：Simpson法求解结果Simps-iont复化梯形公式求解结果匚拈蜘二拦化求槎六、实验中遇到的问题及解决办法在本实验中，遇到的问题主要有两个：(5) 如何将上述的积分算法在MATLAB中实现针对此问题我借鉴了网上以及 课本上的算法的MATLAB实现的程序；(6) 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经

22、验完 成了此项任务。七、实验结论通过以上的测试，我们发现以上算法和程序能够求出积分的比较精确解。八、参考文献1杨大地，王开荣.2006.数值分析.北京:科学出版社何光辉.2008.数值分析实验.数理学院数学实验教学中心3 百度文库，百度知道教师签名年 月 日课程名称数值分析实验项目名 称常微分方程的数值解法实验项目类型验证演示综合设计其他指导教师何光辉成绩V、实验目的：(1) 学会欧拉方法和四阶龙格-库塔方法的使用(2) 设计出相应的算法，编制相应的函数子程序(3) 会用这些函数解决实际问题二、实验内容用欧拉方法和四阶龙格-库塔方法求解微分方程初值问题：y=sin(y)*x，y(0)=10，求

23、y(1)。三、实验原理(1)欧拉法欧拉法是解初值问题的最简单的数值方法。从(9.2)式由于j (%) =y已给定，因而可以算出 v(君)=，(0,j )设X=h充分小，则近似地有：10- w y(X ) = f (X , J ) h00 0记 y. = y(x.) i = O，1，,n从而我们可以取y = y = hf (x，y )1000作为y (x1)的近似值。利用y1及f (x1, yj又可以算出y(x2)的近似值:一般地，在任意点xn+1 = (n + 1)h处y(x)的近似值由下式给出y , = y + hf (x , y )n+1 nn nn+1这就是欧拉法的计算公式，h称为步长。

24、在实际计算时，可将欧拉法与梯形法则相结合，计算公式为:y：+i=yn+hf (气，yn)f (气，yn)+f( xn+i，y;+1hy (k+1)= y+n+1n 2k = 0,1, 2,.(2)四阶龙格-库塔方法 四阶龙格-库塔法求解公式如下:h=儿 + 6(ki 办2 +2k3 + k4)y ,n+1k1 = f (x , y ) n 2 n 2 1) k 3 =fn + n + 2 k 2 / k = f (x + h, y + hk3)k2 = f x + h, y + k四、MATLAB软件实现(1)分别写出欧拉方法和四阶龙格-库塔方法求解微分方程的求解函数%Euler 法，初值 y

25、(a)=c%function y=Euler(f,a,b,c)n=1000;h=(b-a)/n;X=a:h:b;Y=zeros(1,n+1);Y(1)=c;for i=2:n+1x=X(i-1);y=Y(i-1);Y(i)=Y(i-1)+eval(f)*h;end%EAxAuin-eaE-%function y=RK(f,a,b,c)n=1000;h=(b-a)/n;X=a:h:b;Y=zeros(1,n+1);Y(1)=c;for i=1:nx=X(i);y=Y(i);K1=h*eval(f);x=x+h/2;y=y+K1/2;K2=h*eval(f);x=x;y=Y(i)+K2/2;K3=

26、h*eval(f);x=X(i)+h;y=Y(i)+K3;K4=h*eval(f);Y(i+1)=Y(i)+(K1+2火K2+2火K3+K4)/6;end(2)建立界面利用MATLAB中的GUI编程建立如下界面：微分方程数值解(y-f(x3y)请输入函数f 序)请输入初值点xo和所求点X1请输入初值点X。的函数值V(XO)I EljI。法求解结果|口阶RK公式求解结Cnrinht- 沛 圮IT对m时FA对详见程序。五、计算实例、数据、结果、分析下面我们对以上的问题进行测试：输入数据：计算结果如下:六、实验中遇到的问题及解决办法在本实验中，遇到的问题主要有两个：(7) 如何将上述的求解微分方程的

27、算法在MATLAB中实现针对此问题我借鉴了网上以及 课本上的算法的MATLAB实现的程序；(8) 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完 成了此项任务。七、实验结论通过以上的测试，我们发现以上算法和程序能够求出微分方程组的比较精确解。八、参考文献1杨大地，王开荣.2006.数值分析.北京:科学出版社何光辉.2008.数值分析实验.数理学院数学实验教学中心3百度文库，百度知道教师签名年 月 日课程名称数值分析实验项目名 称估计水塔的水流量实验项目类型验证演示综合设计其他指导教师何光辉成绩V一、

28、实验目的：（1）学会对实际问题的分析方法（2）学会利用所学的知识解决实际问题（3）设计出相应的算法，编制相应的应用程序实验内容某居民区，其自来水是有一个圆柱形水塔提供，水塔高12.2m，塔的直径为17.4m，水塔是由水 泵根据水塔中的水位自动加水，一般水泵每天工作两次。按照设计，当水塔中的水位降低至最低水 位，约8.2m时，水泵自动启动加水。当水位升至最高水位，约10.8m时，水泵停止工作。表略。三、实验原理计算中将流量定义为单位时间流出的水的高度乘以水塔横截面积。把时间分成 5段：第1未供水段、水泵开启第1段、第2未供水段、水泵开启第2段、第3未 供水段。先直接对第1、2、3未供水段进行5次

29、曲线拟合。再对得到的曲线分别求 导，取得流速（即单位时间内流出的水的高度）。水泵开启第1、2段，分别在两端 各取两个点，用时刻流速进行拟合得到这两段的流速。流速乘以水塔横截面积就得 到任何时刻的水流量。对其进行分段积分，求和得到一天的总水流量。四、MATLAB软件实现（1）程序function pushbutton1_Callback(hObject, eventdata, handles)% hObject handle to pushbutton1 (see GCBO)% eventdata reserved - to be defined in a future version of M

30、ATLAB% handles structure with handles and user data (see GUIDATA) figure(1);x=0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46636,4995;y=31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.97,35.50,34.45,33.5 0,32.67,31.56,30.81,30.12,29.27,28.42,27.67,26.97,34.75,33.89,33.4

31、0;t=x/3 60 0;%时间单位为小时h=y/3.281;%水位高度为米x1=t(1:10);y1=h(1:10);f1=polyfit(x1,y1,5);t1=0:0.01:t(10);h1=polyval(f1,t1);plot(x1,y1,o,t1,h1,k);xlabel(时间(h);ylabel(水位(m);title(第一阶段供水的时间水位图)function pushbutton2_Callback(hObject, eventdata, handles)% hObject handle to pushbutton2 (see GCBO)% eventdata reserve

32、d - to be defined in a future version of MATLAB% handles structure with handles and user data (see GUIDATA) figure(2);x=0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46636,4995; y=31.75,31.10,30.54,29.94,29.55,28.92,28.50,27.87,27.52,26.97,35.50,34.45,33.5 0,32.67,31.56,30.81,30.

33、12,29.27,28.42,27.67,26.97,34.75,33.89,33.40;t=x/3 60 0;%时间单位为小时h=y/3.281;%水位高度为米x2=t(11:21);y2=h(11:21); f2=polyfit(x2,y2,5); t2=t(11):0.01:t(21);h2=polyval(f2,t2);plot(x2,y2,o,t2,h2,r);xlabel(时间(h);ylabel(水位(m);title(第二阶段供水时段水位图)% Executes on button press in pushbutton4.function pushbutton4_Callba

34、ck(hObject, eventdata,%handles)hObject eventdatahandle to pushbutton4 (see GCBO) reserved - to be defined in a futureversion of MATLABhandlesstructure with handles and user data(see GUIDATA)figure(3);x=0,3316,6635,10619,13937,17921,21240,25223,28543,32284,39435,43318,46636,4995;y=31.75,31.10,30.54,2

35、9.94,29.55,28.92,28.50,27.87,27.52,26.97,35.50,34.45,33.5 0,32.67,31.56,30.81,30.12,29.27,28.42,27.67,26.97,34.75,33.89,33.40;t=x/3600;h=y/3.281;x3=t(22:24);y3=h(22:24);f3=polyfit(x3,y3,5);t3=t(22):0.01:t(24);h3=polyval(f3,t3);plot(x3,y3,o,t3,h3,r);xlabel(时间(h);ylabel(水位(m);title(第三阶段的时间水位图)% Execut

36、es on button press in pushbutton7.function pushbutton7_Callback(hObject, eventdata,%handles)hObject eventdatahandle to pushbutton7 (see GCBO) reserved - to be defined in a futureversion of MATLABhandlesstructure with handles and user data(see GUIDATA)figure(5);x1=0 0.921479 9.308 9.125 8.982 8.814 8

37、.686 8.525 8.388 8.220; 八2.*y1./4;1.843 2.949 3.871 4.978 5.900 7.006 7.928 8.967;y1=9.677 9.v1=pi.*17.4.dv1=-gradient(v1,x1);x2=10.594 12.032 12.954 13.875 14.985 15.903 16.826 17.931 19.037 19.959 20.839;y2=10.820 10.500 10.210 9.936 9.653 9.409 9.180 8.921 8.662 8.433 8.220; v2=pi.*17.4.八2.火y2./4

38、;dv2=-gradient(v2,x2);x3=22.958 23.880 24.986 25.908;y3=10.820 10.591 10.354 10.180;v3=pi.*17.4.八2.火y3./4;dv3=-gradient(v3,x3);x4=x1 x2 x3;v4=dv1 dv2 dv3;x=0:0.001:25.908;v=interp1 (x4,v4,x,spline);plot (x4,v4,+,x,v)(2)建立界面利用MATLAB中的GUI编程建立如下界面： untiledI d I 回 I 详见程序。五、计算实例、数据、结果、分析以上的结果分别如下图：第三未供水时

39、段的时间水位图10.610.5510.510.4S10.4S 10.3510.32 O J5-3520.1间时524.六、实验中遇到的问题及解决办法在本实验中，遇到的问题主要有两个：(9) 如何将上述的建模实现，针对此问题我借鉴了网上以及 课本上的算法的MATLAB实现的程序;(10) 如何将建立界面使得可以随意输入想要求解的相关矩阵后就可以直接求解针对此问题，我通过网上的一些关于MATLAB的GUI设计的相关资料，总结经验完 成了此项任务。七、实验结论通过以上的测试，我们发现以上算法和程序能够求出比较精确的值。八、参考文献1杨大地，王开荣.2006.数值分析.北京:科学出版社何光辉.2008.数值分析实验.数理学院数学实验教学中心3百度文库，百度知道教师签名