数值分析与MATLAB课程 实验程序.docx
- 文档编号:7522509
- 上传时间:2023-01-24
- 格式:DOCX
- 页数:14
- 大小:107.07KB
数值分析与MATLAB课程 实验程序.docx
《数值分析与MATLAB课程 实验程序.docx》由会员分享,可在线阅读,更多相关《数值分析与MATLAB课程 实验程序.docx(14页珍藏版)》请在冰豆网上搜索。
数值分析与MATLAB课程实验程序
第二次实验内容:
首先下载10.22.74.73网站下数值计算课程“实验内容”栏目里的第二章进行学习和实践(其中“二、矩阵的特殊参数运算;三、矩阵的分解;四、矩阵中的一些特殊处理函数”部分暂时不看),然后完成:
1)定义-4pi<=x<=4pi,y1=sinx/x,y2=2sin(x+0.4pi)/(x+0.4pi),在同一坐标系中画出上面两条曲线,两条曲线用不同的颜色,不同的连线标记,不同的顶点标记来加以区分。
(提示:
用plot指令,可用helpplot学习MATLAB的在线帮助)
2)输入如下两个矩阵A和B,对矩阵A和B作关系运算,标识出两矩阵中元素相等的位置,元素值不等的位置,并标识出矩阵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')//作函数y1黑色-
holdon//hold住
plot(x,y2,'red+')//作函数y2红色+
title('曲线')//标注
xlabel('x')
ylabel('y')
holdoff//不用hold了
2.A=[123;-213;-321];
B=[143;328;523];
A==B//这是逻辑运算,符合1,不符合0
A~=B
A<0
3.A=[123;-213;-321];
B=[143;328;523];
A|B
A&B
(B>2)&(B<5)
4.A=[123;-213;-321];
B=[143;328;523];
C=[(A(2,:
))'B(:
3)(A(1,:
))'(A(3,:
))']
运行结果:
1.
2.ans=
101
000
010
ans=
010
111
101
ans=
000
100
100
3.ans=
111
111
111
ans=
111
111
111
ans=
011
100
001
4.C=
-231-3
1822
3331
第二次实验内容:
对于非线性方程
(1)编写M–File函数用二分法求出其在区间[0,0.5]和[2.5,3]内的根,要求函数的最大循环次次为1000次,两个精度均为10-4。
要求打印出最后的根及误差以及循环次数
对于方程
(2)编写M–File函数用迭代法求出其根,初始值为0.5,精度为0.05,最大循环次次为10000次,要求打印出最后的根及误差以及循环次数
1.functionbmz(a,b,epsx,epsy,n)
fork=1:
n
x=(a+b)/2;
if(abs(fun1(x)) fprintf('x=%ff(x)=%f\nxe=%fk=%d\nOK! \n',x,fun1(x),xe,k) return; end iffun1(a)*fun1(x)<0 b=x; else a=x; end xe=abs(b-a)/2; end functiony=fun1(x) y=(x-pi/2)^2-sin(x)-1; 运行结果: bmz(0,0.5,1e-4,1e-4,1000) x=0.394287f(x)=0.000024 xe=0.000244k=11 OK! bmz(2.5,3,1e-4,1e-4,1000) x=2.747314f(x)=0.000053 xe=0.000244k=11 OK! 2.functionitera(x0,eps,n) fork=1: n x=fun2(x0); xe=abs(x-x0); ifxe fprintf('x=%f\nxe=%fk=%d\nOK! \n',x,xe,k) break; else x0=x; end end functiony=fun2(x) y=(x^2-3)/2; 运行结果: itera(0.5,0.05,10000) x=-1.024839 xe=0.049995k=3183 OK! 1、对于下列方程求根 2X1+X2+X3=7 4X1+5X2-X3=11 X1+X2+X3=o (1)用MATLAB语言自己编写一个高斯消去法的程序,用于解上一个方程的根。 (2)选做: 用用MATLAB语言自己编写一个LU分解法的程序,用于解上一个方程的根。 提示: 关于高斯消元法和诺尔当消元法老师上得PPT讲的比较透彻,而且PPT里面也有程序。 functionx=gauss(A,b) n=length(b); a=[Ab]; fork=1: n-1%求上三角矩阵需要注意的是一般写出来的函数是没有常量的。 forj=k+1: n m(j,k)=a(j,k)/a(k,k); fori=1: n+1 a(j,i)=a(j,i)-a(k,i)*m(j,k); end end end x=0;%注意,这一步不能少,少了会得出很离谱的结果 fori=n: -1: 1%回代 x(i)=(a(i,n+1)-sum(a(i,1: n).*x))/a(i,i); end 运行: A=[211;45-1;111]; b=[7;11;0]; >>x=gauss(A,b) x= 7-4-3 对于LU分解,有内部函数lu,我觉得自己编写LU分解不是很重要,但是看看,连连手还是不错的。 有关资料参考老师的PPT--CKCNA03b。 第八页 function[L,U]=myLU(A) %实现对矩阵A的LU分解,L为下三角矩阵 [n,n]=size(A); L=zeros(n,n); U=zeros(n,n); fori=1: n L(i,i)=1; end fork=1: n forj=k: n U(k,j)=A(k,j)-sum(L(k,1: k-1).*U(1: k-1,j)'); end fori=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= 100 210 0.50.166671 U= 211 03-3 001 y= 7 -3 -3 x= 7-4-3 第五次实验内容: 1、此次实验要求每个人都要完成: 编写MATLAB函数用雅可比法解下面方程,要求最大的迭代次数为100,精度控制为||xk+1-xk||∞<0.01,初始向量为零向量。 要求打印出解向量和最终的迭代次数。 2.对于完成实验内容1的同学如果有时间,再编写高斯赛得尔算法完成上面的任务,并且将这两个算法结果作一个比较 对于雅克比,这里就不再写了 2.function[x,k,error]=GaussSaidel(n,A,b,eps,N) fori=1: n x(i)=0; end k=1; whilek<=N error=0; fori=1: n x0=x(i); sum=0; forj=1: n%这里本来要用sum函数,但是在运行的时候老是出现下标有问题,所以改成了for循环 ifj~=i sum=sum+A(i,j)*x(j); end end x(i)=(b(i)-sum)/A(i,i); iferror error=abs(x0-x(i)); end end iferror disp('OK') return; else k=k+1; end end disp('sorry') A=[10-2-1;-210-1;-1-25]; b=[31510]; eps=0.01; N=100; n=3; [x,k,error]=GaussSaidel(n,A,b,eps,N) OK x= 0.99971.99992.9999 k= 5 error= 0.001878 1.任意给定一组数据(n+1)个点;(如下数据仅供参考) 1)、用Lagrange插值法得到n次多项式曲线1; 2)、用最小二乘法进行[n/2]次多项式曲线拟合得到曲线2; 对于第一题,想输出像拟合一样输出一个系数矩阵。 也就是说我最终得到的是一个系数矩阵。 functionp=lagelanri(dx,dy) n=length(dx); symsxyl y=0; fori=1: n l(i)=1; forj=1: n ifj==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); end p=sym2poly(y); 命令窗口: dx=[12345678910]; dy=[3.03.73.94.25.76.67.16.74.53.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,'*') holdon plot(x,y1) plot(x,y2,'red') legend('数据点','插值曲线1','拟合曲线2') holdoff 运行结果: p1= Columns1through8 -0.000100860.0051959-0.114561.4118-10.64950.565-149.56262.32 Columns9through10 -243.3892.4 p2= 0.0037051-0.10150.97566-4.00547.3561-1.24 x 1 2 5 6 7 8 10 13 17 23 f(x) 3.0 3.7 3.9 4.2 5.7 6.6 7.1 6.7 4.5 3.0 自动步长积分求积实验: 把曲线在它的定义域中用自动步长迭代法 我这里的复合梯形是把P175的P1(x)反复的计算。 拉格朗日插值函数和前面一样。 function[sum,k,error]=comptrape(dx,dy,a,b,eps) n=1; h=(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; fork=1: m sum=0; x=a: h: b; y=polyval(p1,x); n=length(x); fori=1: n-1 sum=sum+h*(y(i)+y(i+1))/2; end error=abs(sum1-sum); iferror break; else h=h/2; sum1=sum; end end 命令窗口: a=1;b=23; eps=0.01; dx=[12567810131723]; dy=[3.03.73.94.25.76.67.16.74.53.0]; >>[sum,cishu,e]=comptrape(dx,dy,a,b,eps) sum= 5674.2 cishu= 14 e= 0.0062471 综合练习: 2.任意给定一组数据(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; ift==1 p1=lagelanri(dx,dy); else p1=polyfit(dx,dy,5); end x=a: h: b; y=polyval(p1,x); T1=(b-a)*(y (1)+y (2))/2; fork=1: m sum=0; x=a: h: b; y=polyval(p1,x); n=length(x); fori=1: n-1 sum=sum+h*(y(i)+y(i+1))/2; end error=abs(sum1-sum); iferror ift==1 disp('曲线1OK! ! ! '); else disp('曲线2OK! ! ! '); end return; else h=h/2; sum1=sum; end end 命令窗口的数据如下: clc clear dx=[12345678910]; dy=[3.03.73.94.25.76.67.16.74.53.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,'*') holdon plot(x,y1) plot(x,y2,'red') legend('数据点','插值曲线1','拟合曲线2') holdoff [sum,k,error]=comptr(dx,dy,1,10,0.001,1) [sum,k,error]=comptr(dx,dy,1,10,0.001,2)
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 数值分析与MATLAB课程 实验程序 数值 分析 MATLAB 课程 实验 程序