数值分析与MATLAB课程 实验程序.docx
《数值分析与MATLAB课程 实验程序.docx》由会员分享,可在线阅读,更多相关《数值分析与MATLAB课程 实验程序.docx(14页珍藏版)》请在冰豆网上搜索。
![数值分析与MATLAB课程 实验程序.docx](https://file1.bdocx.com/fileroot1/2023-1/24/3b2dd60c-1065-49c4-b045-4ed3839c74ef/3b2dd60c-1065-49c4-b045-4ed3839c74ef1.gif)
数值分析与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);
ifxefprintf('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);
iferrorerror=abs(x0-x(i));
end
end
iferrordisp('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);
iferrorbreak;
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);
iferrorift==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)