数值分析与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）;
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）
升级会员 