常州大学数值分析作业共六章.docx
- 文档编号:27096049
- 上传时间:2023-06-27
- 格式:DOCX
- 页数:38
- 大小:270.88KB
常州大学数值分析作业共六章.docx
《常州大学数值分析作业共六章.docx》由会员分享,可在线阅读,更多相关《常州大学数值分析作业共六章.docx(38页珍藏版)》请在冰豆网上搜索。
常州大学数值分析作业共六章
第一章:
clear
f=@(x)1/2-x^2/24;
f(0.012)
ans=
0.5000
9.设
,给出计算函数值
的一个合适算法,并在字长m给定的,十进制计算机上给出数值计算结果。
解:
由
得
10.字长为5的十进制计算机上计算
和
,并与
的精确值
1.79比较,说明差异存在理由,其中
,
。
解:
字长为5时的误差很大,这是因为设置的字长有限,就不可避免的使舍入误差不断积累。
把字长改为9时,误差已经大幅度减小。
这说明,加大字长可以显著减小误差。
11.举例介绍数组矩阵常见运算。
解:
举例如下
12.对任意给定的实数a、b、c、试编写Matlab程序,求方程
的根。
clear
a=input('输入a=');
b=input('输入b=');
c=input('输入c=');
d=b^2-4*a*c
ifb>=0
x1=(-b-d^0.5)/2*a
x2=-2*c/(d^0.5+b)
elseb<0
x1=-2*c/(d^0.5-b)
x2=(-b+d^0.5)/2*a
end
输入a=2
输入b=2
输入c=1
d=
-4
x1=
-2.0000-2.0000i
x2=
-0.5000+0.5000i
输入a=1
输入b=2
输入c=1
d=
0
x1=
-1
x2=
-1
输入a=1
输入b=5
输入c=2
d=
17
x1=
-4.5616
x2=
-0.4384
解:
利用教材例11的方法:
当b>0时,
,
。
当b<0时,
,
。
13.利用
及
,给出一个计算π的方法,根据此方法编写程序,给出π的至少有10位有效数字的近似值。
clear
x=3^(-0.5);
y=atan(x);%精确值%
s=0;%计算值%
fork=1:
2:
100;
s=s+(-1)^((k+1)/2)*(x^k)/k;
err=y-s;
m=abs(err)
ifm<=1e-11
break
end
end
pi=6*s
解:
根据题中所给公式,容易得到:
14.分别利用下式给出计算ln2的近似方法,编写相应程序并比较算法运行情况。
function[x,iternum,flag]=jacobi(A,b,x0,delta,max1)
%检验输入参数,初始化
ifnargin<2,error('moreaugmentsareneeded');end
ifnargin<3,x0=zeros(size(b));end
ifnargin<4,delta=1e-13;end
ifnargin<5,max1=100;end
ifnargin>5,error('incorrectnumberofinput');end
n=length(b);x=0*b;flag=0;iternum=0;
%用Jacobi迭代法解方程组
fork=1:
max1
iternum=iternum+1;
fori=1:
n
ifabs(A(i,i)) error('A(i,i)equaltozero,dividedbyzero'); end x(i)=(b(i)-A(i,[1: i-1,i+1: n])*x0([1: i-1,i+1: n]))/A(i,i); end err=norm(x-x0); relerr=err/(norm(x)+eps); x0=x; if(err flag=1; break; end end ifflag==1 disp('TheJacobimethodconverges.') x=x; else disp(['TheJacobimethoddoesnotconvergewith'... num2str(max1),'iterations']) end return clear x=1/3; s=0; fork=1: 2: 100; s=s+((x)^k)/k; end y=2*s%计算值% g=log (2)%真实值% err=y-g%绝对误差% y= 0.6931 g= 0.6931 err= -2.2204e-16 clear x=1; s=0; fork=1: 100; s=s+((-1)^(k+1))*(x^k)/k; end y=s%计算值% g=log (2)%真实值% err=y-g%绝对误差% y= 0.6882 g= 0.6931 err= -0.0050 解: 由运行结果可知, 方法二的绝对误差比方法一的误差要小得多。 这是因为方法一给出的计算公式含有相近数相减项,损失了有效数字。 而方法二给出的计算公式避免了相近数相减,具有较好的精度。 第二章: 20. (1)用Jacobi迭代法解方程组AX=b. function[x,iternum,flag]=gseid(A,b,x0,delta,max1) %检验输入参数,初始化 ifnargin<2,error('moreaugmentsareneeded');end ifnargin<3,x0=zeros(size(b));end ifnargin<4,delta=1e-13;end ifnargin<5,max1=100;end ifnargin>5,error('incorrectnumberofinput');end n=length(b);flag=0;iternum=0; %用Gauss-Seidel迭代法解方程组 fork=1: max1 iternum=iternum+1; x=x0; fori=1: n ifabs(A(i,i)) error('A(i,i)equaltozero,dividedbyzero'); end x0(i)=(b(i)-A(i,[1: i-1,i+1: n])*x0([1: i-1,i+1: n]))/A(i,i); end err=norm(x-x0); relerr=err/(norm(x0)+eps); if(err flag=1; break; end end ifflag==1 disp('TheGauss-seidelmethodconverges.') x=x0; else disp(['TheGauss-seidelmethoddoesnotconvergewith'... num2str(max1),'iterations']) end return A=[12-2;111;221];b=[1;1;1]; [x,iternum,flag]=jacobi(A,b) 输出结果为: TheJacobimethodconverges. ans=-3 3 1 iternum=4 flag=1 20. (1)用Gauss-Seidel迭代法解方程组AX=b. A=[12-2;111;221];b=[1;1;1];[x,iternum,flag]=gseid(A,b) 输出结果为: TheGauss-seidelmethoddoesnotconvergewith100iterations ans=1.0e+31* -9.1905 9.2222 -0.0634 iternum=4 flag=1 A= 12-2 111 221 b= 1 1 1 function[L,U]=Crout(A) %检验输入参数,初始化 b=size(A);n=b (1); ifb (1)~=b (2) error('MatrixAshouldbeaSquarematrix.'); end ifn~=rank(A) error('MatrixAshouldbeFULLRANK.'); end L=zeros(n,n);U=zeros(n,n); fori=1: n U(i,i)=1; end %将矩阵A进行Crout分解 fork=1: n fori=k: n temp_sum=0; fort=1: k-1 temp_sum=temp_sum+L(i,t)*U(t,k); end L(i,k)=A(i,k)-temp_sum; end forj=k+1: n temp_sum=0; fort=1: k-1 temp_sum=temp_sum+L(k,t)*U(t,j); end U(k,j)=(A(k,j)-temp_sum)/L(k,k); end end return 3. (1)将矩阵A进行Crout分解 A=[1020;0111;20-11;0011] A= 1020 0111 20-11 0011 [L,U]=Crout(A) L= 1.0000000 01.000000 2.00000-5.00000 001.00001.2000 U= 1.000002.00000 01.00001.00001.0000 001.0000-0.2000 0001.0000 function[L,U]=Doolittle(A) %检验输入参数,初始化 b=size(A);n=b (1); ifb (1)~=b (2) error('MatrixAshouldbeaSquarematrix.'); end ifn~=rank(A) error('MatrixAshouldbeFULLRANK.'); end L=zeros(n,n);U=zeros(n,n); fori=1: n L(i,i)=1; end %将矩阵A进行Doolittle分解 fork=1: n forj=k: n temp_sum=0; fort=1: k-1 temp_sum=temp_sum+L(k,t)*U(t,j); end U(k,j)=A(k,j)-temp_sum; end fori=k+1: n temp_sum=0; fort=1: k-1 temp_sum=temp_sum+L(i,t)*U(t,k); end L(i,k)=(A(i,k)-temp_sum)/U(k,k); end end return 3. (2)将矩阵A进行Doolittle分解。 A=[1020;0111;20-11;0011] A= 1020 0111 20-11 0011 [L,U]=Doolittle(A) L= 1.0000000 01.000000 2.000001.00000 00-0.20001.0000 U= 1.000002.00000 01.00001.00001.0000 00-5.00001.0000 0001.2000 function[x]=improvecholesky(A,b,n) %检验输入参数,初始化 L=zeros(n,n);D=diag(n,0);S=L*D; fori=1: n L(i,i)=1; end fori=1: n forj=1: n if(eig(A)<=0)|(A(i,j)~=A(j,i)) disp('MatrixAshouldbeaPositive-definitematrix.'); break; end end end D(1,1)=A(1,1); %用改进平方根法解方程组 fori=2: n forj=1: i-1 S(i,j)=A(i,j)-sum(S(i,1: j-1)*L(j,1: j-1)'); L(i,1: i-1)=S(i,1: i-1)/D(1: i-1,1: i-1); end D(i,i)=A(i,i)-sum(S(i,1: i-1)*L(i,1: i-1)'); end y=zeros(n,1);x=zeros(n,1); fori=1: n y(i)=(b(i)-sum(L(i,1: i-1)*D(1: i-1,1: i-1)*y(1: i-1)))/D(i,i); end fori=n: -1: 1 x(i)=y(i)-sum(L(i+1: n,i)'*x(i+1: n)); end return 7.用改进平方根法解方程组AX=b function[L,U,x]=pursue(a,b,c,f) n=length(b);u (1)=b (1); fori=2: n if(u(i-1)~=0) l(i-1)=a(i-1)/u(i-1); u(i)=b(i)-l(i-1)*c(i-1); else break; end end L=eye(n)+diag(l,-1); U=diag(u)+diag(c,1); x=zeros(n,1);y=x;y (1)=f (1); fori=2: n y(i)=f(i)-l(i-1)*y(i-1); end if(u(n)~=0) x(n)=y(n)/u(n); end fori=n-1: -1: 1 x(i)=(y(i)-c(i)*x(i+1))/u(i); end return A=[41-10;13-10;-1-152;0024]b=[1;0;0;0] A= 41-10 13-10 -1-152 0024 b= 1 0 0 0 n=4 [x]=improvecholesky(A,b,n) n= 4 x= 0.2821 -0.0769 0.0513 -0.0256 8. (2)用追赶法求解方程组. a=[1-1-1-1];b=[22222];c=[-1-1-1-1];f=[10000]; [L,U,x]=pursue(a,b,c,f) L= 1.00000000 0.50001.0000000 0-0.40001.000000 00-0.62501.00000 000-0.72731.0000 U= 2.0000-1.0000000 02.5000-1.000000 001.6000-1.00000 0001.3750-1.0000 00001.2727 x= 0.3571 -0.2857 -0.2143 -0.1429 第三章: 1.设节点x0=0,x1=π/8,x2=π/4,x3=3π/8,x4=π/2,试适当选取上述节点,用Lagrange插值法分别构造cosx在区间[0,π/2]上的一次、二次、四次差值多项式P1(x),P2(x)和P4(x),并分别计算P1(π/3),P2(π/3)和P4(π/3)。 function[A,Y]=lagrange(x,y,x0) %检验输入参数 ifnargin<2||nargin>3 error('IncorrectNumberofInputs'); end iflength(x)~=length(y) error('Thelengthofxmustbeequaltoitofy'); end m=length(x);n=m-1;L=zeros(m,m); %计算基本插值多项式的系数 fori=1: n+1 C=1; forj=1: n+1 ifi~=j ifabs(x(i)-x(j)) error('therearetwosamenodes'); end C=conv(C,poly(x(j)))/(x(i)-x(j)); end end L(i,: )=C; end %计算Lagrange插值多项式的系数 A=y*L; %计算f(x0)的近似值 ifnargin==3 Y=polyval(A,x0); end A=fliplr(A); return x=[0,pi/4,pi/2]; y=cos(x); x0=pi/3; [A,Y]=lagrange(x,y,x0); P2=vpa(poly2sym(A),3) Y P2= x^2-0.109*x-0.336 Y= 0.5174 x=[0,pi/8,pi/4,3*pi/8,pi/2]; y=cos(x); x0=pi/3; [A,Y]=lagrange(x,y,x0); P4=vpa(poly2sym(A),3) Y P4= x^4+0.00282*x^3-0.514*x^2+0.0232*x+0.0287 Y= 0.5001 x=[pi/8,3*pi/8]; y=cos(x); x0=pi/3; [A,Y]=lagrange(x,y,x0); P1=vpa(poly2sym(A),3) Y P1= 1.19*x-0.689 Y= 0.4729 7.根据列表函数,选取适当的节点,用逐次线性插值法给出三次差值多项式在2.8处的值。 x 0.0 1.0 2.0 3.0 4.0 f(x) 0.50 1.25 2.75 3.50 2.75 function[T,y0]=aitken(x,y,x0,T0) ifnargin==3 T0=[]; end n0=size(T0,1); m=max(size(x)); n=n0+m; T=zeros(n,n+1); T(1: n0,1: n0+1)=T0; T(n0+1: n,1)=x; T(n0+1: n,2)=y; ifn0==0 i0=2; else i0=n0+1; end fori=i0: n forj=3: i+1 T(i,j)=fun(T(j-2,1),T(i,1),T(j-2,j-1),T(i,j-1),x0); end end y0=T(n,n+1); return function[y]=fun(x1,x2,y1,y2,x) y=y1+(y2-y1)*(x-x1)/(x2-x1); return %选取0、1、3、4四个节点,求三次插值多项式 x=[0,1,3,4]; y=[0.5,1.25,3.5,2.75]; x0=2.8; [T,y0]=aitken(x,y,x0) T= 00.5000000 1.00001.25002.600000 3.00003.50003.29993.23000 4.00002.75002.07502.28503.4190 y0= 3.4190 16. function[C,D,Y]=newpoly(x0,y0,x) %检验输入参数 ifnargin<2|nargin>3 error('IncorrectNumberofInputs'); end iflength(x0)~=length(y0) error('Thelengthofx0mustbeequaltoitofy0'); end n=length(x0);D=zeros(n,n);D(: 1)=y0'; %计算差商表 forj=2: nfork=j: n ifabs(x0(k)-x0(k-j+1)) error('DividedbyZero,therearetwonodesarethesame'); end D(k,j)=(D(k,j-1)-D(k-1,j-1))/(x0(k)-x0(k-j+1)); end end %计算Newton插值多项式的系数 C=D(n,n); fork=(n-1): -1: 1 C=conv(C,poly(x0(k))); m=length(C); C(m)=C(m)+D(k,k); end ifnargin==3 Y=polyval(C,x); end C=fliplr(C); return 选取适当的函数y=f(x)和插值节点,编写Matlab程序,分别给出利用Lagrange插值方法、Newton插值方法确定的插值多项式,并将函数y=f(x)的插值多项式及插值余项的图形画在同一坐标系中,观测节点变化对插值余项的影响。 第四章: 6.已知列表函数,用最小二乘法求形如 的拟合函数。 x 1.0 2.0 3.0 4.0 5.0 f(x) 1.222 2.984 5.466 8.902 13.592 x=[12345]; y=[1.2222.9845.4668.90213.592]; [a,b]=ec(x,y) 输出结果为: a= 1.3205 b= 0.2786 plot(x,y,'b*',x,1*x.*exp(0.2*x),'r-') function[a,b]=ec(x,y) Y=log(y)'; A=zeros(5,3); fori=1: 5 A(i,1)=1; A(i,2)=log(x(i)); A(i,3)=i; end c=inv(A'*A)*(A'*Y); a=exp(c (1)); b=c(3); fori=1: 5 y=a*x.*exp(b*x); end return 7.已知人体表面积S和身高h,体重w有近似关系 ,试根据身高、体重及相应的人体表面积的一组观测值(hi,wi,Si)(i=0,1,2,,n)来估计参数a0a1a2的大小。 function[a0,a1,a2]=ec2(h,w) S=log(s)';N=length(h); A=z
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 常州 大学 数值 分析 作业 共六章
![提示](https://static.bdocx.com/images/bang_tan.gif)