数值分析作业

1、MATLAB数学实验报告 姓名： 劳杰林 学号： 指导老师： 易昆南 专业班级： 统计1002 一解线性方程组的直接法学号 班级统计1002姓名 劳杰林 指导教师易昆南实验题目Gauss消去法、列主元消去法、LU分解法评 分1、设计（实习）目的：1. 学习gauss消去原理及过程。2.利用matlab直接求解线性方程组。2、实验内容：利用matlab，编写gauss消去法、列主元消去法及LU分解法的程序，直接解线性方程组。3.详细设计： gauss消去法： function x=guass(A,b)n=rank(A);for k=1:n-1 for i=k+1:n for j=k+1:n if

2、 abs(A(k,k)1e-6 warning(分母为零不能计算!); return; else A(i,j)=A(i,j)-(A(i,k)/A(k,k)*A(k,j); end end b(i)=b(i)-(A(i,k)/A(k,k)*b(k); A(i,k)=0; endendformat long ex(n)=b(n)/A(n,n);for k=n-1:-1:1 sum=0; for j=k+1:n sum=sum+A(k,j)*x(j); end x(k)=(b(k)-sum)/A(k,k);end列主元消去法：function x=gaussxiaoqumethod(A,b)n=le

3、ngth(b);x=zeros(n,1);c=zeros(1,n);t=0;for i=1:n-1 max=abs(A(i,i); m=i; for j=i+1:n if max x=guass(2.5 2.3 -5.1;5.3 9.6 1.5;8.1 1.7 -4.3,3.7;3.8;5.5)x =4.9263e-001 2.3071e-001 -4.5055e-001列主元消去法： x=gaussxiaoqumethod(0.101 2.301 3.555;-1.347 3.712 4.623;-2.835 1.072 5.643,1.183;2.137;3.035)x = -3.7065

4、e-001 1.7579e-0023.4636e-001LU分解法： l u=lufj(1 2 3;2 5 2;3 1 5)l = 1 0 0 2 1 0 3 -5 1u = 1 2 3 0 1 -4 0 0 -245：实验总结1.用matlab求解线性方程组问题有着很强的实用性，它能帮助我们简化繁琐又难以计算的数学问题。2.在设计过程中我发现先前看似无味的数学也有几丝美妙，这让我更加坚定了学习这门课程的信心与动力，相信通过深入地学习，今后定能有所收获。3.设计过程中，我深刻体会到与同学交流是十分重要的。比如说有时候一个小小的错误自己总是检查不出来，但让同学帮忙看一下马上就查出来了。 姓名：劳

5、杰林 2012年3月21日二牛顿插值学号 班级统计1002姓名 劳杰林 指导教师易昆南实验题目牛顿插值、复化梯形公式、复化辛普森公式、节点加密复化梯形公式评 分1、设计（实习）目的：1.利用matlab直接求解积分。2、实验内容：1.牛顿插值法2.复化梯形公式3.复化辛普森公式4.节点加密复化梯形公式3.详细设计： 牛顿插值法： clear allx0=-2 -1 0 1 3;y0=-56 -16 -2 -2 4;n=length(x0);A=zeros(n);for i=1:n A(i,1)=y0(i);endfor j=2:n for k=j:n A(k,j)=(A(k,j-1)-A(k-

6、1,j-1)/(x0(k)-x0(k-j+1); end end a=diag(A); syms x f; s=0; p=1; for i=1:(n-1) for j=1:i p=p*(x-x0(j); end s=s+a(i+1)*p; p=1; end f=a(1)+s 复化梯形公式：format longa=input(a=);b=input(b=);n=input(n=);h=(b-a)/n;x0=a:h:b;syms x;y=input(y=);y0=subs(y,x0);s1=0;for i=2:n-1 s1=s1+y0(i);ends=h/2*(y0(1)+y0(n)+2*s1)

7、;disp(利用复化梯形公式求解积分数值为,num2str(s)复化辛普森公式：a=input(a=);b=input(b=);n=input(n=);h=(b-a)/(n);x0=a:h/2:b;syms x;y=input(y=);y0=subs(y,x0);s2=0;s3=0;for i=2:2:2*n-1 s2=s2+y0(i);endfor i=3:2:2*n-1 s3=s3+y0(i);endt=h/6*(y0(1)+y0(n)+2*s3+4*s2);disp(利用复化Simpson公式求得的积分数值为,num2str(t)节点加密复化梯形公式：format longa=input

8、(a=);b=input(b=);n=input(n=);h=(b-a)/(2*n);x0=a:h:b;syms x;y=input(y=);y0=subs(y,x0);s1=0; s2=0; for i=3:2:(2*n-2) s1=s1+y0(i);endfor j=2:2:(2*n-1) s2=s2+y0(j);ends=h/2*(y0(1)+y0(2*n)+2*s1)+h*s2disp(利用节点加密复化梯形公式求解积分数值为,num2str(s)4：实验结果：牛顿插值法： f = 40*x - 13*(x + 1)*(x + 2) + 2*x*(x + 1)*(x + 2) + 24复

9、化梯形公式：a=0b=1n=50y=x利用复化梯形公式求解积分数值为0.4802复化辛普森公式：a=0b=1n=50y=x利用复化Simpson公式求得的积分数值为0.4851节点加密复化梯形公式：a=0b=1n=25y=xs = 0.0000利用节点加密复化梯形公式求解积分数值为0.461 5：实验总结1.用matlab求解线性方程组问题有着很强的实用性，它能帮助我们简化繁琐又难以计算的数学问题。2.在设计过程中我发现先前看似无味的数学也有几丝美妙，这让我更加坚定了学习这门课程的信心与动力，相信通过深入地学习，今后定能有所收获。3.设计过程中，我深刻体会到与同学交流是十分重要的。比如说有时候

10、一个小小的错误自己总是检查不出来，但让同学帮忙看一下马上就查出来了。 姓名：劳杰林 2012年4月5日三、Romberg法求积分学号 班级统计1002姓名 劳杰林 指导教师易昆南实验题目编程求积分，dx，01sinxxdx评 分1、设计（实习）目的：1 了解MATLAB在实际问题中的应用2 通过实践加深对这门语言中M文件的了解3 掌握运用matlab解决实际数学问题2、实验内容：1）用复化梯形公式计算2）用romberg求积方法计算积分dx3)用romberg算法计算积分I=3.详细设计： 1.） function y=f(x)y=e-x2;function T_n=F_H_T(a,b,n)h

11、=(b-a)/n;x=a:h:b;for k=0:n x(k+1)=a+k*h; if x(k+1)=0 x(k+1)=10(-10); endendT_1=h/2*(f(x(1)+f(x(n+1);for i=2:n F(i)=h*f(x(i);endT_2=sum(F);T_n=T_1+T_2;Romberg算法：function s,n=rombg2(a,b,eps1)s=10;s0=0;k=3;t=;t(1,1)=(b-a)*(f(b)+f(a)/2;while abs(s-s0)eps1 h=(b-a)/2(k-1); w=0; if h=0 for i=1:2(k-1)-1 w=w

12、+f(a+i*h); end t(k,1)=h*(f(a)/2+f(b)/2+w)/2; for l=2:k for i=1:(k-l+1) t(i,l)=(4(l-1)*t(i+1,l-1)-t(i,l-1)/(4(l-1)-1); end end s=t(1,k); s0=t(1,k-1); k=k+1; n=k; else s=s0; n=-k; endend2） function y=f(x)y=(2*exp(-x)/sqrt(pi);3） function y=f(x)y=(sin(x)/x;4：实验结果：1.） T_n=F_H_T(0,1,10)T_n = 0.74622） s,n

13、=rombg2(0,1,5e-6)s = 0.3566n = 83.） s,n=rombg2(1e-8,1,5e-8)s = 0.4730n = 85：实验总结1.用matlab求解线性方程组问题有着很强的实用性，它能帮助我们简化繁琐又难以计算的数学问题。2.在设计过程中我发现先前看似无味的数学也有几丝美妙，这让我更加坚定了学习这门课程的信心与动力，相信通过深入地学习，今后定能有所收获。3.设计过程中，我深刻体会到与同学交流是十分重要的。比如说有时候一个小小的错误自己总是检查不出来，但让同学帮忙看一下马上就查出来了。 姓名：劳杰林 2012年4月02日四、常微分求解学号 班级统计1002姓名

14、劳杰林 指导教师易昆南实验题目编写四个解线性方程组的程序评 分1、设计（实习）目的：4 了解MATLAB在实际问题中的应用5 通过实践加深对这门语言中M文件的了解6 掌握运用matlab解决实际数学问题2、实验内容：用欧拉方法，休恩方法，泰勒方法，龙格库塔方法0 3解常微分方程，初值 ,y(0)=1,h=1,1/2,1/4,1/8,并画出图形。3.详细设计： 1.）欧拉法functionxout,yout=euler(ypfun,xspan,y0,h)f=ypfun;,x=xspan(1):h:xspan(2);n=length(x);y=zeros(1,n);y(1)=y0;for i=1:

15、n-1 y(i+1)=y(i)+h*feval(f,x(i),y(i);endxout=x;yout=y; function f=myfun(x,y)f=(x-y)/2;2.）休恩法 functionxout,yout=eulerpro(ypfun,xspan,y0,h)f=ypfun;x=xspan(1):h:xspan(2);n=length(x);y=zeros(1,n);y(1)=y0;for i=1:n-1 yp=y(i)+h*feval(f,x(i),y(i); yc=y(i)+h*feval(f,x(i+1),yp); y(i+1)=(yp+yc)/2;endxout=x;you

16、t=y;function f=myfun(x,y)f=(x-y)/2;3）泰勒方法function x y=taylor(h) x0=0;y0=1;ya0=(x0-y0)/2;yb0=(2-x0+y0)/4;yc0=(-2+x0-y0)/8;yd0=(2-x0+y0)/16;x(1)=h;y(1)=y0+h*ya0+(h2)/2*yb0+(h3)/6*yc0+(h4)/24*yd0;ya(1)=(x(1)-y(1)/2;yb(1)=(2-x(1)+y(1)/4;yc(1)=(-2+x(1)-y(1)/8;yd(1)=(2-x(1)+y(1)/16;for n=1:(3/h)-1 x(n+1)=

17、h*(n+1); y(n+1)=y(n)+h*ya(n)+(h2)/2*yb(n)+(h3)/6*yc(n)+(h4)/24*yd(n); ya(n+1)=(x(n+1)-y(n+1)/2;yb(n+1)=(2-x(n+1)+y(n+1)/4;yc(n+1)=(-2+x(n+1)-y(n+1)/8;yd(n+1)=(2-x(n+1)+y(n+1)/16;end4）龙格库塔法function x y =lungkuta(h)x0=0;y0=1;k10=h*(x0-y0)/2;k20=h*(x0+0.5-y0-k10/2)/2;k30=h*(x0+0.5-y0-k20/2)/2;k40=h*(x0

18、+0.5-y0-k30/2)/2;y(1)=y0+h*(k10+2*k20+2*k30+k40)/6;x(1)=h;k1(1)=h*(x(1)-y(1)/2;k2(1)=h*(x(1)+0.5-y(1)-k1(1)/2)/2;k3(1)=h*(x(1)+0.5-y(1)-k2(1)/2)/2;k4(1)=h*(x(1)+0.5-y(1)-k3(1)/2)/2;y(2)=y(1)+h*(k1(1)+2*k2(1)+2*k3(1)+k4(1)/6;for n=1:(3/h)-1 x(n+1)=h*(n+1); y(n+1)=h*y(n)+(k1(n)+2*k2(n)+2*k3(n)+k4(n)/6

19、; k1(n+1)=h*(x(n+1)-y(n+1)/2; k2(n+1)=h*(x(n+1)+0.5-y(n+1)-k1(n+1)/2)/2; k3(n+1)=h*(x(n+1)+0.5-y(n+1)-k2(n+1)/2)/2; k4(n+1)=h*(x(n+1)+0.5-y(n+1)-k3(n+1)/2)/2;end4：实验结果：1.） x1,y1=euler(myfun,0 3,1,1)x2,y2=euler(myfun,0 3,1,1/2)x3,y3=euler(myfun,0 3,1,1/4)x4,y4=euler(myfun,0 3,1,1/8)plot(x1,y1,-,x2,y2

20、,r,x3,y3,b,x4,y4,g)x = 0 1 2 3y = 1.0000 0.5000 0.7500 1.3750x = 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000y = 1.0000 0.7500 0.6875 0.7656 0.9492 1.2119 1.5339x = Columns 1 through 12 0 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500 Column 13 3.0000y = Columns 1 throug

21、h 12 1.0000 0.8750 0.7969 0.7598 0.7585 0.7887 0.8464 0.9281 1.0308 1.1520 1.2892 1.4406 Column 13 1.6043x = Columns 1 through 12 0 0.1250 0.2500 0.3750 0.5000 0.6250 0.7500 0.8750 1.0000 1.1250 1.2500 1.3750 Columns 13 through 24 1.5000 1.6250 1.7500 1.8750 2.0000 2.1250 2.2500 2.3750 2.5000 2.6250

22、 2.7500 2.8750 Column 25 3.0000y = Columns 1 through 12 1.0000 0.9375 0.8867 0.8469 0.8174 0.7976 0.7868 0.7845 0.7902 0.8033 0.8234 0.8500 Columns 13 through 24 0.8829 0.9214 0.9654 1.0144 1.0682 1.1265 1.1889 1.2552 1.3252 1.3986 1.4753 1.5549 Column 25 1.63742.） x1,y1=eulerpro(myfun,0 3,1,1)x2,y2

23、=eulerpro(myfun,0 3,1,1/2)x3,y3=eulerpro(myfun,0 3,1,1/4)x4,y4=eulerpro(myfun,0 3,1,1/8)plot(x1,y1,-,x2,y2,r,x3,y3,b,x4,y4,g)x = 0 1 2 3y = 1.0000 0.8750 1.1719 1.7324x = 0 0.5000 1.0000 1.5000 2.0000 2.5000 3.0000y = 1.0000 0.8438 0.8311 0.9305 1.1176 1.3731 1.6821x = Columns 1 through 12 0 0.2500

24、0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500 Column 13 3.0000y = Columns 1 through 12 1.0000 0.8984 0.8381 0.8141 0.8222 0.8587 0.9201 1.0037 1.1068 1.2271 1.3626 1.5115 Column 13 1.6723x = Columns 1 through 12 0 0.1250 0.2500 0.3750 0.5000 0.6250 0.7500 0.8750 1.0000 1.1250

25、 1.2500 1.3750 Columns 13 through 24 1.5000 1.6250 1.7500 1.8750 2.0000 2.1250 2.2500 2.3750 2.5000 2.6250 2.7500 2.8750 Column 25 3.0000y = Columns 1 through 12 1.0000 0.9434 0.8977 0.8624 0.8368 0.8203 0.8124 0.8125 0.8202 0.8350 0.8565 0.8842 Columns 13 through 24 0.9178 0.9570 1.0013 1.0506 1.10

26、44 1.1625 1.2247 1.2907 1.3602 1.4332 1.5092 1.5883 Column 251.67013）： x1 y1=taylor(1)x2 y2=taylor(1/2)x3 y3=taylor(1/4)x4 y4=taylor(1/8)plot(x1,y1,-,x2,y2,r,x3,y3,g,x4,y4,b)x1 = 1 2 3y1 = 0.0000 1.2084 1.3738x2 = 0.0000 1.0000 1.0000 2.0000 2.0000 3.0000y2 = 0.0000 0.6558 0.5024 1.0246 1.6256 1.915

27、9x3 = Columns 1 through 6 0.0000 0.0000 0.0000 1.0000 1.0000 1.0000 Columns 7 through 12 1.0000 2.0000 2.0000 2.0000 2.0000 3.0000y3 = Columns 1 through 6 0.8125 0.7229 0.3750 0.0794 0.3871 0.8097 Columns 7 through 12 1.4770 1.5855 1.5184 1.0557 1.9650 1.7015x4 = Columns 1 through 6 0.0000 0.0000 0.

28、0000 0.0000 0.0000 0.0000 Columns 7 through 12 0.0000 1.0000 1.0000 1.0000 1.0000 1.0000 Columns 13 through 18 1.0000 1.0000 1.0000 2.0000 2.0000 2.0000 Columns 19 through 24 2.0000 2.0000 2.0000 2.0000 2.0000 3.0000y4 = Columns 1 through 6 0.6133 0.6014 0.0613 0.0323 0.1696 0.0840 Columns 7 through

29、 12 0.3207 0.4163 0.6442 0.1253 0.1410 0.9040 Columns 13 through 18 0.8933 1.2965 1.2707 1.7712 1.7307 1.5297 Columns 19 through 24 1.0865 1.6254 1.3503 1.9607 1.2746 1.2613 4) x1 y1 =lungkuta(1) x2 y2 =lungkuta(1/2) x3 y3 =lungkuta(1/4) x4 y4 =lungkuta(1/8) plot(x1,y1,-,x2,y2,r,x3,y3,g,x4,y4,b) x1 = 1 2 3y1 = 0.0000 1.2188 1.2228x2 = 0.0000 1.0000 1.0000 2.0000 2.0000 3.0000y2 = 0.7500 0.3609 0.4637 0.1974 0.5500 0.7