常州大学数值分析作业第三章.docx
- 文档编号:4653116
- 上传时间:2023-05-07
- 格式:DOCX
- 页数:42
- 大小:110.29KB
常州大学数值分析作业第三章.docx
《常州大学数值分析作业第三章.docx》由会员分享,可在线阅读,更多相关《常州大学数值分析作业第三章.docx(42页珍藏版)》请在冰点文库上搜索。
常州大学数值分析作业第三章
第一章:
clear
f=@(x)1/2-x^2/24;
f(0.012)
ans=
0.5000
9.设
,给出计算函数值
的一个合适算法,并在字长m给定的,十进制计算机上给出数值计算结果。
解:
由
得
10.字长为5的十进制计算机上计算
和
,并与
的精确值
1.0075376410479比较,说明差异存在理由,其中
,
。
clear
f=@(x)digit(digit(exp(x)-1,5)/x,5);
g=@(x)digit(digit(1,5)+digit(x/2,5)+digit...
(digit(x^2,5)/6,5)+digit(digit(x^3,5)/24,5),5);
exc=1.0075376410479;
f(0.015)
g(0.015)
err1=f(0.015)-exc
err2=g(0.015)-exc
ans=
1.0075
ans=
1.0075
err1=
-3.7641e-05
err2=
-3.7641e-05
%%编写m文件使用digit函数设置字长%%
functiony=digit(x,m)
k=max(size(x));
y=x;
fori=1:
k
ifx(i)<0
sign=-1;
else
sign=1;
end
x(i)=abs(x(i));
p=0;
ifx(i)<0.1&x(i)>eps
whilex(i)<0.1
x(i)=x(i)*10;
p=p-1;
end
end
ifx(i)>=1
whilex(i)>=1
x(i)=x(i)/10;
p=p+1;
end
end
y(i)=round(x(i)*10^m)/10^m;
y(i)=sign*y(i)*10^p;
end
return
f=@(x)digit(digit(exp(x)-1,9)/x,9);
g=@(x)digit(digit(1,9)+digit(x/2,9)+digit...
(digit(x^2,9)/6,9)+digit(digit(x^3,9)/24,9),9);
err1=f(0.015)-exc
err2=g(0.015)-exc
err1=
-1.0479e-09
err2=
-1.0479e-09
解:
字长为5时的误差很大,这是因为设置的字长有限,就不可避免的使舍入误差不断积累。
把字长改为9时,误差已经大幅度减小。
这说明,加大字长可以显著减小误差。
11.举例介绍数组矩阵常见运算。
解:
举例如下
A.*B
ans=
1234
10121416
27303336
52566064
A*B
ans=
30303030
70707070
110110110110
150150150150
clear
A=[1:
4;5:
8;9:
12;13:
16]
B=[1,1,1,1;2,2,2,2;3,3,3,3;4,4,4,4]
A’
A=
1234
5678
9101112
13141516
B=
1111
2222
3333
4444
ans=
15913
261014
371115
481216
A^2
ans=
90100110120
202228254280
314356398440
426484542600
A.^2
ans=
14916
25364964
81100121144
169196225256
A-B
ans=
0123
3456
6789
9101112
A+B
ans=
2345
78910
12131415
17181920
12.对任意给定的实数a、b、c、试编写Matlab程序,求方程
的根。
解:
利用教材例11的方法:
当b>0时,
,
。
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
当b<0时,
,
。
输入a=2
输入b=2
输入c=1
d=
-4
x1=
-2.0000-2.0000i
x2=
-0.5000+0.5000i
输入a=1
输入b=5
输入c=2
d=
17
x1=
-4.5616
x2=
-0.4384
输入a=1
输入b=2
输入c=1
d=
0
x1=
-1
x2=
-1
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的近似方法,编写相应程序并比较算法运行情况。
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
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
解:
由运行结果可知,
方法二的绝对误差比方法一的误差要小得多。
这是因为方法一给出的计算公式含有相近数相减项,损失了有效数字。
而方法二给出的计算公式避免了相近数相减,具有较好的精度。
20.分别用Jacobi迭代法、Gauss-seidel迭代法解方程组
解:
Jacobi迭代法收敛
,Gauss-seidel迭代法不收敛。
27.编写LU分解法、改进平方根法、追赶法的Matlab程序,并进行相关数值实验。
3.将矩阵
进行Doolittle和Crout分解
解:
Doolittle分解:
结果如下,程序见后面。
Crout分解:
结果如下,程序见后面。
7.用改进平方根法解方程组
解:
结果如下,程序见后面。
8
(2).用追赶法求解方程组
解:
结果如下,程序见后面。
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 A=[12-2;111;221] b=[1;1;1] A= 12-2 111 221 b= 1 1 1 jacobi(A,b) TheJacobimethodconverges. ans= -3 3 1 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 gseid(A,b) TheGauss-seidelmethoddoesnotconvergewith100iterations ans= 1.0e+31* -9.1905 9.2222 -0.0634 A=[12-2;111;221] b=[1;1;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 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]=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进行Doolittle分解 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 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 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 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 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=[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 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 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]=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-i),x0); end end return function[y]=fun(x1,x2,y1,y2,x) y=y1+(y2-y1)*(x-x1)/(x2-x1); return x=[01]; y=[0.51.25]; x0=2
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 常州 大学 数值 分析 作业 第三