常州大学数值分析作业第三章Word文件下载.docx
- 文档编号:6819763
- 上传时间:2023-05-07
- 格式:DOCX
- 页数:42
- 大小:110.29KB
常州大学数值分析作业第三章Word文件下载.docx
《常州大学数值分析作业第三章Word文件下载.docx》由会员分享,可在线阅读,更多相关《常州大学数值分析作业第三章Word文件下载.docx(42页珍藏版)》请在冰点文库上搜索。
p=p-1;
ifx(i)>
=1
whilex(i)>
x(i)=x(i)/10;
p=p+1;
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);
-1.0479e-09
字长为5时的误差很大,这是因为设置的字长有限,就不可避免的使舍入误差不断积累。
把字长改为9时,误差已经大幅度减小。
这说明,加大字长可以显著减小误差。
11.举例介绍数组矩阵常见运算。
举例如下
A.*B
1234
10121416
27303336
52566064
A*B
30303030
70707070
110110110110
150150150150
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=
5678
9101112
13141516
B=
1111
2222
3333
4444
15913
261014
371115
481216
A^2
90100110120
202228254280
314356398440
426484542600
A.^2
14916
25364964
81100121144
169196225256
A-B
0123
3456
6789
A+B
2345
78910
12131415
17181920
12.对任意给定的实数a、b、c、试编写Matlab程序,求方程
的根。
利用教材例11的方法:
当b>
0时,
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<
x1=-2*c/(d^0.5-b)
x2=(-b+d^0.5)/2*a
当b<
输入a=2
输入b=2
输入c=1
d=
-4
x1=
-2.0000-2.0000i
x2=
-0.5000+0.5000i
输入a=1
输入b=5
输入c=2
17
-4.5616
-0.4384
0
-1
13.利用
及
,给出一个计算π的方法,根据此方法编写程序,给出π的至少有10位有效数字的近似值。
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
pi=6*s
根据题中所给公式,容易得到:
14.分别利用下式给出计算ln2的近似方法,编写相应程序并比较算法运行情况。
x=1;
s=0;
s=s+((-1)^(k+1))*(x^k)/k;
y=s%计算值%
g=log
(2)%真实值%
err=y-g%绝对误差%
y=
0.6882
g=
0.6931
err=
-0.0050
x=1/3;
s=s+((x)^k)/k;
y=2*s%计算值%
-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'
3,x0=zeros(size(b));
4,delta=1e-13;
5,max1=100;
ifnargin>
5,error('
incorrectnumberofinput'
n=length(b);
x=0*b;
flag=0;
iternum=0;
%用Jacobi迭代法解方程组
max1
iternum=iternum+1;
fori=1:
n
ifabs(A(i,i))<
error('
A(i,i)equaltozero,dividedbyzero'
x(i)=(b(i)-A(i,[1:
i-1,i+1:
n])*x0([1:
n]))/A(i,i);
err=norm(x-x0);
relerr=err/(norm(x)+eps);
x0=x;
if(err<
delta)||(relerr<
delta)
flag=1;
break;
ifflag==1
disp('
TheJacobimethodconverges.'
)
x=x;
else
disp(['
TheJacobimethoddoesnotconvergewith'
...
num2str(max1),'
iterations'
])
A=[12-2;
111;
221]
b=[1;
1;
1]
12-2
111
221
b=
1
jacobi(A,b)
TheJacobimethodconverges.
-3
3
function[x,iternum,flag]=gseid(A,b,x0,delta,max1)
%用Gauss-Seidel迭代法解方程组
x=x0;
x0(i)=(b(i)-A(i,[1:
relerr=err/(norm(x0)+eps);
TheGauss-seidelmethodconverges.'
TheGauss-seidelmethoddoesnotconvergewith'
gseid(A,b)
TheGauss-seidelmethoddoesnotconvergewith100iterations
1.0e+31*
-9.1905
9.2222
-0.0634
function[L,U]=Crout(A)
b=size(A);
n=b
(1);
ifb
(1)~=b
(2)
MatrixAshouldbeaSquarematrix.'
ifn~=rank(A)
MatrixAshouldbeFULLRANK.'
end
L=zeros(n,n);
U=zeros(n,n);
U(i,i)=1;
%将矩阵A进行Crout分解
fori=k:
n
temp_sum=0;
fort=1:
k-1
temp_sum=temp_sum+L(i,t)*U(t,k);
L(i,k)=A(i,k)-temp_sum;
forj=k+1:
temp_sum=temp_sum+L(k,t)*U(t,j);
U(k,j)=(A(k,j)-temp_sum)/L(k,k);
end
A=[1020;
0111;
20-11;
0011]
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
%将矩阵A进行Doolittle分解
[L,U]=Doolittle(A)
2.000001.00000
00-0.20001.0000
00-5.00001.0000
0001.2000
function[x]=improvecholesky(A,b,n)
D=diag(n,0);
S=L*D;
L(i,i)=1;
forj=1:
if(eig(A)<
=0)|(A(i,j)~=A(j,i))
MatrixAshouldbeaPositive-definitematrix.'
D(1,1)=A(1,1);
%用改进平方根法解方程组
fori=2:
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);
D(i,i)=A(i,i)-sum(S(i,1:
i-1)*L(i,1:
i-1)'
y=zeros(n,1);
x=zeros(n,1);
y(i)=(b(i)-sum(L(i,1:
i-1)*D(1:
i-1)*y(1:
i-1)))/D(i,i);
fori=n:
-1:
1
x(i)=y(i)-sum(L(i+1:
n,i)'
*x(i+1:
n));
A=[41-10;
13-10;
-1-152;
0024]
0;
0]
41-10
13-10
-1-152
0024
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)
u
(1)=b
(1);
if(u(i-1)~=0)
l(i-1)=a(i-1)/u(i-1);
u(i)=b(i)-l(i-1)*c(i-1);
L=eye(n)+diag(l,-1);
U=diag(u)+diag(c,1);
y
(1)=f
(1);
y(i)=f(i)-l(i-1)*y(i-1);
if(u(n)~=0)
x(n)=y(n)/u(n);
fori=n-1:
1
x(i)=(y(i)-c(i)*x(i+1))/u(i);
a=[1-1-1-1];
b=[22222];
c=[-1-1-1-1];
f=[10000];
[L,U,x]=pursue(a,b,c,f)
1.00000000
0.50001.0000000
0-0.40001.000000
00-0.62501.00000
000-0.72731.0000
2.0000-1.0000000
02.5000-1.000000
001.6000-1.00000
0001.3750-1.0000
00001.2727
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>
IncorrectNumberofInputs'
iflength(x)~=length(y)
Thelengthofxmustbeequaltoitofy'
m=length(x);
n=m-1;
L=zeros(m,m);
%计算基本插值多项式的系数
n+1
C=1;
ifi~=j
ifabs(x(i)-x(j))<
therearetwosamenodes'
C=conv(C,poly(x(j)))/(x(i)-x(j));
L(i,:
)=C;
%计算Lagrange插值多项式的系数
A=y*L;
%计算f(x0)的近似值
ifnargin==3
Y=polyval(A,x0);
A=fliplr(A);
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);
P2=vpa(poly2sym(A),3)
Y
P2=
x^2-0.109*x-0.336
0.5174
x=[0,pi/8,pi/4,3*pi/8,pi/2];
[A,Y]=lagrange(x,y,x0);
P4=vpa(poly2sym(A),3)
P4=
x^4+0.00282*x^3-0.514*x^2+0.0232*x+0.0287
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
function[T]=aitken(x,y,x0,T0)
ifnargin==3
T0=[];
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;
n,2)=y;
ifn0==0
i0=2;
i0=n0+1;
fori=i0:
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);
function[y]=fun(x1,x2,y1,y2,x)
y=y1+(y2-y1)*(x-x1)/(x2-x1);
x=[01];
y=[0.51.25];
x0=2
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 常州 大学 数值 分析 作业 第三