数值分析与MATLAB课程 实验程序.docx

1、数值分析与MATLAB课程 实验程序第二次实验内容：首先下载10.22.74.73网站下数值计算课程“实验内容”栏目里的第二章进行学习和实践（其中“二、矩阵的特殊参数运算；三、矩阵的分解；四、矩阵中的一些特殊处理函数”部分暂时不看），然后完成：1） 定义-4pi=x=4pi，y1=sinx/x, y2=2sin(x+0.4pi)/(x+0.4pi)，在同一坐标系中画出上面两条曲线，两条曲线用不同的颜色，不同的连线标记，不同的顶点标记来加以区分。（提示：用plot指令，可用help plot学习MATLAB的在线帮助） 2） 输入如下两个矩阵 A 和 B，对矩阵 A 和 B 作关系运算，标识出两

2、矩阵中元素相等的位置，元素值不等的位置，并标识出矩阵 A 中所有小于 0 的元素。3) 对上面的矩阵 A 和 B 作逻辑“ 或”、“ 与”运算，并标识出矩阵 B 中所有大于 2 并小于 5 的元素位置。4) 利用矩阵的寻址方式，取出上面矩阵A中第二行元素作为新矩阵的第一列，B矩阵的第3列元素作为为新矩阵的第二列元素，取出矩阵A的第一行，第三行元素作为新矩阵的第三列和第四列元素， 1. x=-4*pi:0.1:4*pi;y1=sin(x)./x;y2=2*sin(x+0.4*pi)./(x+0.4*pi);/数组所以是用点除。应该是数组和数组的运算用点plot(x,y1,-black)/作函数y

3、1黑色-hold on/hold住plot(x,y2,red+) /作函数y2红色+title(曲线)/标注xlabel(x)ylabel(y)hold off/不用hold了2. A=1 2 3;-2 1 3;-3 2 1;B=1 4 3;3 2 8;5 2 3;A=B/这是逻辑运算，符合1，不符合0A=BA2)&(B5)4.A=1 2 3;-2 1 3;-3 2 1;B=1 4 3;3 2 8;5 2 3;C=(A(2,:) B(:,3) (A(1,:) (A(3,:)运行结果：1. 2. ans = 1 0 1 0 0 0 0 1 0ans = 0 1 0 1 1 1 1 0 1ans

4、= 0 0 0 1 0 0 1 0 03. ans = 1 1 1 1 1 1 1 1 1ans = 1 1 1 1 1 1 1 1 1ans = 0 1 1 1 0 0 0 0 14. C = -2 3 1 -3 1 8 2 2 3 3 3 1第二次实验内容：对于非线性方程（1） 编写M File函数用二分法求出其在区间0，0.5和2.5，3内的根，要求函数的最大循环次次为1000次，两个精度均为10-4。要求打印出最后的根及误差以及循环次数对于方程（2） 编写M File函数用迭代法求出其根，初始值为0.5，精度为0.05，最大循环次次为10000次，要求打印出最后的根及误差以及循环次数1

5、.function bmz(a,b,epsx,epsy,n)for k=1:n x=(a+b)/2; if (abs(fun1(x)epsy )|(abs(b-a)/2)epsx) fprintf(x=%f f(x)=%f nxe=%f k=%dnOK!n,x,fun1(x),xe,k) return; end if fun1(a)*fun1(x)0 b=x; else a=x; end xe=abs(b-a)/2;end function y=fun1(x)y=(x-pi/2)2-sin(x)-1;运行结果：bmz(0,0.5,1e-4,1e-4,1000)x=0.394287 f(x)=0

6、.000024 xe=0.000244 k=11OK!bmz(2.5,3,1e-4,1e-4,1000)x=2.747314 f(x)=0.000053 xe=0.000244 k=11OK!2. function itera(x0,eps,n)for k=1:n x=fun2(x0); xe=abs(x-x0); if xe x=gauss(A,b)x = 7 -4 -3对于LU分解，有内部函数lu，我觉得自己编写LU分解不是很重要，但是看看，连连手还是不错的。有关资料参考老师的PPT-CKCNA03b。第八页function L,U=myLU(A)%实现对矩阵A的LU分解，L为下三角矩阵n

7、,n=size(A);L=zeros(n,n);U=zeros(n,n);for i=1:n L(i,i)=1;endfor k=1:n for j=k:n U(k,j)=A(k,j)-sum(L(k,1:k-1).*U(1:k-1,j); end for i=k+1:n L(i,k)=(A(i,k)-sum(L(i,1:k-1).*U(1:k-1,k)/U(k,k); end end运行结果；L = 1 0 0 2 1 0 0.5 0.16667 1U = 2 1 1 0 3 -3 0 0 1y = 7 -3 -3x = 7 -4 -3第五次实验内容：1、 此次实验要求每个人都要完成：编写M

8、ATLAB 函数用雅可比法解下面方程，要求最大的迭代次数为100，精度控制为 |xk+1-xk|0.01，初始向量为零向量。 要求打印出解向量和最终的迭代次数。2. 对于完成实验内容1 的同学如果有时间，再编写高斯赛得尔算法完成上面的任务，并且将这两个算法结果作一个比较对于雅克比，这里就不再写了2. function x,k,error=GaussSaidel(n,A,b,eps,N)for i=1:n x(i)=0;endk=1;while k=N error=0; for i=1:n x0=x(i); sum=0; for j=1:n%这里本来要用sum函数，但是在运行的时候老是出现下标有

9、问题，所以改成了for循环 if j=i sum=sum+A(i,j)*x(j); end end x(i)=(b(i)-sum)/A(i,i); if errorabs(x0-x(i) error=abs(x0-x(i); end end if erroreps disp(OK) return; else k=k+1; endenddisp(sorry)A=10 -2 -1;-2 10 -1;-1 -2 5;b=3 15 10;eps=0.01;N=100;n=3;x,k,error=GaussSaidel(n,A,b,eps,N)OKx = 0.9997 1.9999 2.9999k =

10、5error = 0.0018781. 任意给定一组数据（n+1）个点；(如下数据仅供参考)1）、用Lagrange插值法得到n次多项式曲线1；2）、用最小二乘法进行n/2次多项式曲线拟合得到曲线2；对于第一题，想输出像拟合一样输出一个系数矩阵。也就是说我最终得到的是一个系数矩阵。function p=lagelanri(dx,dy)n=length(dx);syms x y ly=0;for i=1:n l(i)=1; for j=1:n if j=i else l(i)=l(i)*(x-dx(j)/(dx(i)-dx(j); end end l(i)=l(i)*dy(i); y=y+l(i

11、);endp=sym2poly(y);命令窗口：dx=1 2 3 4 5 6 7 8 9 10;dy=3.0 3.7 3.9 4.2 5.7 6.6 7.1 6.7 4.5 3.0;p1=lagelanri(dx,dy)x=1:0.01:10;y1=polyval(p1,x);p2=polyfit(dx,dy,5)y2=polyval(p2,x);plot(dx,dy,*)hold onplot(x,y1)plot(x,y2,red)legend(数据点,插值曲线1,拟合曲线2)hold off运行结果：p1 = Columns 1 through 8 -0.00010086 0.005195

12、9 -0.11456 1.4118 -10.649 50.565 -149.56 262.32 Columns 9 through 10 -243.38 92.4p2 =0.0037051 -0.1015 0.97566 -4.0054 7.3561 -1.24x12567810131723f(x)3.03.73.94.25.76.67.16.74.53.0自动步长积分求积实验：把曲线在它的定义域中用自动步长迭代法我这里的复合梯形是把P175的P1（x）反复的计算。拉格朗日插值函数和前面一样。function sum,k,error=comptrape(dx,dy,a,b,eps)n=1;h=

13、(b-a)/n;sum1=0;m=100;p1=lagelanri(dx,dy);x=a:h:b;y=polyval(p1,x);T1=(b-a)*(y(1)+y(2)/2;for k=1:m sum=0; x=a:h:b; y=polyval(p1,x); n=length(x); for i=1:n-1 sum=sum+h*(y(i)+y(i+1)/2; end error=abs(sum1-sum); if error sum,cishu,e=comptrape(dx,dy,a,b,eps)sum = 5674.2cishu = 14e = 0.0062471综合练习：2. 任意给定一组

14、数据（n+1）个点；(如下数据仅供参考)2.用Lagrange插值法得到n次多项式曲线1；3.用最小二乘法进行n/2次多项式曲线拟合得到曲线2；4.把2条曲线与原数据一起画在同一坐标下进行比较；5.把2条曲线在它的定义域中用自动步长迭代法进行积分。比较计算的结果。function sum,k,error=comptr(dx,dy,a,b,eps,t)n=1;h=(b-a)/n;sum1=0;m=100;if t=1 p1=lagelanri(dx,dy);else p1=polyfit(dx,dy,5);endx=a:h:b;y=polyval(p1,x);T1=(b-a)*(y(1)+y(2

15、)/2;for k=1:m sum=0; x=a:h:b; y=polyval(p1,x); n=length(x); for i=1:n-1 sum=sum+h*(y(i)+y(i+1)/2; end error=abs(sum1-sum); if erroreps if t=1 disp(曲线1OK!); else disp(曲线2OK!); end return; else h=h/2; sum1=sum; endend命令窗口的数据如下：clccleardx=1 2 3 4 5 6 7 8 9 10;dy=3.0 3.7 3.9 4.2 5.7 6.6 7.1 6.7 4.5 3.0;p1=lagelanri(dx,dy);p2=polyfit(dx,dy,5);x=1:0.01:10;y1=polyval(p1,x);y2=polyval(p2,x);plot(dx,dy,*)hold onplot(x,y1)plot(x,y2,red)legend(数据点,插值曲线1,拟合曲线2 )hold offsum,k,error=comptr(dx,dy,1,10,0.001,1)sum,k,error=comptr(dx,dy,1,10,0.001,2)