上海大学课程实验报告数值代数.docx
- 文档编号:30049095
- 上传时间:2023-08-04
- 格式:DOCX
- 页数:40
- 大小:1.47MB
上海大学课程实验报告数值代数.docx
《上海大学课程实验报告数值代数.docx》由会员分享,可在线阅读,更多相关《上海大学课程实验报告数值代数.docx(40页珍藏版)》请在冰豆网上搜索。
上海大学课程实验报告数值代数
课程实验报告
COURSEPAPER
课程名称:
数值代数与计算方法
课程号:
08305114
授课教师:
学号:
姓名:
所属:
计算机科学与工程
打印时间:
2015
评语:
题目一:
算法:
1、对于问题一:
2、对于问题二:
直接编写递归函数程序,算出三个差分方程的10个近似值
程序:
1)
%main.m
clc;
clearall;
a=1;
b=-2;
c=-3;
[x1,x2]=roots(a,b,c)
%roots.m
function[x1,x2]=roots(a,b,c)
d=sqrt(b*b-4*a*c);
ifd>0
x1=(-2*c)/(b+d);
x2=(-b-d)/(2*a);
elseifd<0
x1=(-b+d)/(2*a);
x2=(-2*c)/(b-d);
end
2)
%solu.m
function[X,R,P,Q]=solu(X0,R0,P0,P1,Q0,Q1)
X
(1)=X0;
R
(1)=R0;
P
(1)=P0;
P
(2)=P1;
Q
(1)=Q0;
Q
(2)=Q1;
fori=1:
9
X(i+1)=X(i)/2;
X(i)=X(i+1);
end
fori=1:
9
R(i+1)=R(i)/2;
R(i)=R(i+1);
end
fori=3:
10
P(i)=3/2*P(i-1)-1/2*P(i-2);
P(i-1)=P(i);
P(i-2)=P(i-1);
end
fori=3:
10
Q(i)=5/2*Q(i-1)-Q(i-2);
Q(i-1)=Q(i);
Q(i-2)=Q(i-1);
end
%Xn.m
clc;
clearall;
X0=1;
R0=0.994;
P0=1;
P1=0.497;
Q0=1;
Q1=0.497;
[X,R,P,Q]=solu(X0,R0,P0,P1,Q0,Q1);
fori=1:
10
x(i)=i;
end
plot(x,X-R,'r*');
holdon
结果:
1)
x1=
3
x2=
-1
2)
题目二:
算法:
1、
此题我们以第一个式子
为例,求出
g(Xn+1)=(2/(2+3*Xn-Xn.^3))^(1/2)
2、利用二分法:
3、利用牛顿迭代法,并根据题目提示的立方根算法:
此题我们以第三个式子为例:
,设A=71/3,求出x的3次方便是近似值。
程序:
1、
%fixpt.m
function[k,p,err,P]=fixpt(p0,tol,max1)
P
(1)=p0;
fork=2:
max1
P(k)=g(P(k-1));
err=abs(P(k)-P(k-1));
relerr=err/(abs(P(k))+eps);
p=P(k);
if(err break; end end ifk==max1 disp('maximumnumberofiterationsexceeded') end P=P'; %g.m functiony=g(x) y=sqrt(2/(2+3*x-x.^3)); %iteration.m clc; clearall; p0=0.5; tol=1.0e-9; max1=20; [k,p,err,P]=fixpt(p0,tol,max1); 2、 二分法: %bisect.m function[c,err,yc,k,an,cn,bn,ycn]=bisect(a,b,delta) ya=f(a); yb=f(b); an (1)=a; bn (1)=b; ifya*yb>0 disp('NoRootinthisdistrict! '); end max1=1+round((log(b-a)-log(delta))/log (2)); fork=1: max1 c=(a+b)/2; cn(k)=c; yc=f(c); ycn(k)=yc; ifyc==0 a=c; an(k+1)=a; b=c; bn(k+1)=b; elseifyb*yc>0 b=c; bn(k+1)=b; yb=yc; else a=c; an(k+1)=a; ya=yc; end ifb-a break; end end c=(a+b)/2; cn(max1)=c; err=abs(b-a); yc=f(c); ycn(max1)=yc; fori=1: 9 k(i)=i-1; end %f.m functionh=f(x) h=x*sin(x)-1; %main.m clc; clearall; a=0; b=2; delta=0.000001; [c,err,yc,k,an,cn,bn,ycn]=bisect(a,b,delta); 4、牛顿迭代法求解: %newton.m function[p0,err,k,y]=newton(A,p0,delta,epsilon,max1) fork=1: max1 p1=(2.*p0+A/p0.^2)/3; err=abs(p1-p0); relerr=2*err/(abs(p1)+delta); p0=p1; y=f(p0,A); if(err end %f.m functiony=f(x,A) y=x.^3-A; %main.m clc; clearall; A=(7).^(1/3); p0=2; delta=0.001; epsilon=0.000001; max1=30; [p0,err,k,y]=newton(A,p0,delta,epsilon,max1); x=p0^3; vpa(x,10) 结果: 1、 2、二分法 3、 即71/3的近似值为1.912931183 题目三: 算法: ▲回代算法: 程序: %backsub.m functionX=backsub(A,B) n=length(B); X=zeros(n,1); X(n)=B(n)/A(n,n); fork=n-1: -1: 1 X(k)=(B(k)-A(k,k+1: n)*X(k+1: n))/A(k,k); end %main.m clearall; clc; A=[ 3,-2,1,-1; 0,4,-1,2; 0,0,2,3; 0,0,0,5 ]; B=[8,-3,11,15]'; X=backsub(A,B) det(A) 结果: X= 2 -2 1 3 ans= 120 即x1=2,x2=-2,x3=1,x4=3; |A|=120. 题目四: 求解方程组: 算法: 先将方程组系数矩阵进行上三角变换,再进行回代,最终求出方程组解。 上三角变换的算法: 高斯消元法,AX=B,设增广矩阵为[A|B],进行初等行变换,通过选取主元行并消元,构成使增广矩阵的A矩阵构成上三角矩阵。 选主元策略常用偏序选主元策略: 首先检查位于主对角线或主对角线下方第p列的所有元素,确定行k,它的元素的绝对值最大,即 |akp|=max{|app|,|ap+1p|,…,|aN-1p|,|aNp|} 如果k>p,则交换第k行和第p行。 程序: %uptrbk.m functionX=uptrbk(A,B) [NN1]=size(A); X=zeros(N,1); C=zeros(1,N+1); Aug=[AB]; forp=1: N-1 [Y,j]=max(abs(Aug(p: N,p))); C=Aug(p,: ); Aug(p,: )=Aug(j+p-1,: ); Aug(j+p-1,: )=C; ifAug(p,p)==0 'Awassingular.NOuniquesolution' break end fork=p+1: N m=Aug(k,p)/Aug(p,p); Aug(k,p: N+1)=Aug(k,p: N+1)-m*Aug(p,p: N+1); end end X=backsub(Aug(1: N,1: N),Aug(1: N,N+1)); %main.m clearall; clc; A=[2,4,-6; 1,5,3; 1,3,2]; B=[-4,10,5]'; X=uptrbk(A,B) det(A) 结果: X= -3 2 1 ans= 18 即x1=-3,x2=2,x3=1;|A|=18. 题目五: 算法: 如下线性方程组: 1、雅可比迭代法 2、高斯-赛德尔迭代法: 在此题中我们分别用三组初始值,P1=[x1=0,x2=0,x3=0];P2=[x1=1,x2=2,x3=3];P3=[x1=5,x2=5,x3=5]; 分别用雅可比和高斯进行迭代,得出结果并比较。 程序: %jacobi.m//雅可比迭代函数 functionX=jacobi(A,B,P,delta,max1) N=length(B); fork=1: max1 forj=1: N X(j)=(B(j)-A(j,[1: j-1,j+1: N])*P([1: j-1,j+1: N]))/A(j,j); end err=abs(norm(X'-P)); relerr=err/(norm(X)+eps); P=X'; if(err break end end X=X'; %gseid.m//高斯-赛德尔迭代函数 functionX=gseid(A,B,P,delta,max1) N=length(B); fork=1: max1 forj=1: N ifj==1 X (1)=(B (1)-A(1,2: N)*P(2: N))/A(1,1); elseifj==N X(N)=(B(N)-A(N,1: N-1)*(X(1: N-1))')/A(N,N); else X(j)=(B(j)-A(j,1: j-1)*X(1: j-1)'-A(j,j+1: N)*P(j+1: N))/A(j,j); end end err=abs(norm(X'-P)); relerr=err/(norm(X)+eps); P=X'; if(err break end end X=X'; %main.m clearall; clc; delta=0.0000001; max1=100; A=[ 1,0,1; -1,1,0; 1,2,-3 ]; B=[2,0,0]'; P1=[0,0,0]'; P2=[1,2,3]'; P3=[5,5,5]'; X1=jacobi(A,B,P1,delta,max1); X1=vpa(X1',10) X2=jacobi(A,B,P2,delta,max1); X2=vpa(X2',10) X3=jacobi(A,B,P3,delta,max1); X3=vpa(X3',10) X4=gseid(A,B,P1,delta,max1); X4=vpa(X4',10) X5=gseid(A,B,P2,delta,max1); X5=vpa(X5',10) X6=gseid(A,B,P3,delta,max1); X6=vpa(X6',10) 结果: X1= [1.001606818,1.004722712,1.003011523] X2= [1.0022603,0.9955518438,0.9943430273] X3= [0.9935727281,0.9811091525,0.9879539061] X4= [0,0,0] X5= [3.0,3.0,3.0] X6= [5.0,5.0,5.0] 正确结果为: x1=x2=x3=1,发现雅可比迭代对于此方程组有效。 题目六: 程序: %sin.m clearall; clc; x=-1: 0.1: 1; P_5=x-(x.^3)/factorial(3)+(x.^5)/factorial(5); P_7=x-(x.^3)/factorial(3)+(x.^5)/factorial(5)-(x.^7)/factorial(7); P_9=x-(x.^3)/factorial(3)+(x.^5)/factorial(5)-(x.^7)/factorial(7)+(x.^9)/factorial(9); f=sin(x); plot(x,f,x,P_5,x,P_7,x,P_9); gridon; holdon; x=-1: 0.2: 1; P_5=x-(x.^3)/factorial(3)+(x.^5)/factorial(5); P_7=x-(x.^3)/factorial(3)+(x.^5)/factorial(5)-(x.^7)/factorial(7); P_9=x-(x.^3)/factorial(3)+(x.^5)/factorial(5)-(x.^7)/factorial(7)+(x.^9)/factorial(9); f=sin(x); D=horzcat(x',P_5',P_7',P_9',f') 结果: 1、 2.xP5(x)P7(x)P9(x)sin(x) 题目七: 此题我们以课本P155页为例 算法: 程序: %Rd.m function[B,D,I]=Rd(P,x) N=length(P); B(N)=P(N); fork=N-1: -1: 1 B(k)=P(k)+B(k+1)*x; end D(N-1)=(N-1)*P(N); fork=N-1: -1: 2 D(k-1)=(k-1)*P(k)+D(k)*x; end I(N+1)=P(N)/(N); fork=N: -1: 2 I(k)=P(k-1)/(k-1)+I(k+1)*x; end I (1)=I (2)*x; %main.m clearall; clc; P=[1.28,-0.4,0.2,-0.02]; x=4; [B,D,I]=Rd(P,x); disp(B (1)); disp(D (1)); disp(I (1)); 结果: 经对比,结果与题目所给答案一致。 题目八: 算法: 第一问我们用拉格朗日多项式对题目给出的数据进行插值拟合;第二问利用题目七的算法,将第一问求出来的拉格朗日多项式进行x=1到x=6的积分再除5求得平均值;第三问用plot函数画出函数图像。 拉格朗日插值多项式: 程序: %lagran.m function[C,L]=lagran(X,Y) w=length(X); n=w-1; L=zeros(w,w); fork=1: n+1 V=1; forj=1: n+1 ifk~=j V=conv(V,poly(X(j)))/(X(k)-X(j)); end end L(k,: )=V; end C=Y*L %main.m clearall; clc; X=[1,2,3,4,5,6]; Y=[66,66,65,64,63,63]; x=6; [C,L]=lagran(X,Y); H=fliplr(C); [I]=Rd(H,x); aver=(I (1)-I(6))/6 x=1: 0.1: 6; y=C (1)*x.^5+C (2)*x.^4+C(3)*x.^3+C(4)*x.^2+C(5)*x.^1+C(6); plot(x,y,'r-'); holdon; gridon; 结果: 得P(x)=0.0167+(-0.2917)x+2.0000*x2+(-6.7083)x3+9.9833*x4+61.0000*x5; 平均值为64.4569. 题目九: 算法: 最小二乘法拟合曲线: 程序: %lsline.m function[A,B]=lsline(X,Y) xmean=mean(X); ymean=mean(Y); sumx2=(X-xmean)*(X-xmean)'; sumxy=(Y-ymean)*(X-xmean)'; A=sumxy/sumx2; B=ymean-A*xmean; %main.m clc; clear; X=0.2: 0.2: 1.0; Y1=[3.6,7.3,10.9,14.5,18.2]; Y2=[5.3,10.6,15.9,21.2,26.4]; [A1,B1]=lsline(X,Y1) [A2,B2]=lsline(X,Y2) 结果: A1= 18.2000 B1= -0.0200 A2= 26.4000 B2= 0.0400 即F1k=18.2000*xk-0.0200;F2k=26.4000*xk+0.0400. 题目十: 算法: 我们以(a)、(c)题为例,利用精度为O(h2)中心差分公式: 程序: %difflim.m function[L,n]=difflim(f,x,max,th) h=th; H (1)=h; D (1)=(feval(f,x+h)-feval(f,x-h))/(2*h); E (1)=0; forn=1: 2 h=h/10; H(n+1)=h; D(n+1)=(feval(f,x+h)-feval(f,x-h))/(2*h); E(n+1)=abs(D(n+1)-D(n)); R(n+1)=2*E(n+1)-(abs(D(n+1))+abs(D(n))+eps); end n=2; while(E(n)>E(n+1))&&(n h=h/10; H(n+2)=h; D(n+2)=(feval(f,x+h)-feval(f,x-h))/(2*h); E(n+2)=abs(D(n+2)-D(n+1)); R(n+2)=2*E(n+2)-(abs(D(n+2))+abs(D(n+1))+eps); n=n+1; end n=length(D)-1; L=[H'D'E']; %main.m clearall; clc; max=40; th=1; x1=1/sqrt(3); f1=@(x)60*x^45-32*x^33+233*x^5-47*x^2-77; x2=1/sqrt (2); f2=@(x)sin(cos(1/3)); [L1,n]=difflim(f1,x1,max,th); L1=vpa(L1,14) [L2,n]=difflim(f2,x2,max,th); L2=vpa(L2,14) 结果: L1=[h][Dk][E=|Dk+1-Dk|]L2=[h][Dk][E] [1.0,24149615818.548,0][1.0,0,0] [0.1,82.963075326944,24149615735.585][0.1,0,0] [0.01,75.251162467466,7.7119128594781][0.01,0,0] [0.001,75.174271349987,0.076891117479505] [0.0001,75.173502461752,0.00076888823485888] [0.00001,75.173494772685,0.0000076890671607543] [0.000001,75.173494700209,0.000000072475373258385] [0.0000001,75.173494664682,0.000000035527136788005] [0.00000001,75.173495162062,0.00000049737991503207] 题目十一: 以下列式子为例,分别用两种方式求积分: 算法: 程序: %traprl.m functions=traprl(f,a,b,m) h=(b-a)/m; s=0; fork=1: (m-1) x=a+h*k; s=s+feval(f,x); end s=h*(feval(f,a)+feval(f,b))/2+h*s; %simprl.m functions=simprl(f,a,b,m) h=(b-a)/(2*m); s1=0; s2=0; fork=1: (m-1) x=a+2*k*h; s1=s1+feval(f,x); end fork=1: m x=a+2*(k-1)*h; s2=s2+feval(f,x); end s=h*(feval(f,a)+feval(f,b)+2*s1+4*s2)/3; %main.m clearall; clc; f=@(x)(1+x.^2).^(-1); a=-1; b=1; s_t=traprl(f,a,b,18257); s_s=simprl(f,a,b,256); S_t=vpa(s_t,10) S_v=vpa(s_s,10) 结果: S1_t= 1.570796326 S1_v= 1.570791241 即组合梯形公式积分结果为: 1.570796326 组合辛普森公式积分结果为: 1.570791241.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 上海 大学 课程 实验 报告 数值 代数
![提示](https://static.bdocx.com/images/bang_tan.gif)