数值分析实验Word文档下载推荐.docx
《数值分析实验Word文档下载推荐.docx》由会员分享,可在线阅读,更多相关《数值分析实验Word文档下载推荐.docx(29页珍藏版)》请在冰豆网上搜索。
代码如下:
sum0=0;
forn=1:
50000
y=(-1)^(n+1)/(2*n-1);
sum1=y+sum1;
pi_1=4*sum0;
pi_2=4*sum1;
error=pi_2-pi_1;
sum0=sum1;
%d,pi_2=%0.8f,Theerror:
n,pi_2,error);
ifabs(error)<
10^(-4)
求得结果是:
3.画图y=sinx与其泰勒展开
x=0:
pi/100:
2*pi;
y=sin(x);
y1=0;
y2=0;
y3=0;
fori=0:
2
y1=y1+(-1)^(i)*x.^(2*i+1)/factorial(2*i+1);
5
y2=y2+(-1)^(i)*x.^(2*i+1)/factorial(2*i+1);
10
y3=y3+(-1)^(i)*x.^(2*i+1)/factorial(2*i+1);
plot(x,y,'
*r'
x,y1,'
b'
x,y2,'
-g'
x,y3,'
k'
)
axis([02*pi-1.51.5])
看出n=2发散,n=10几乎与y=sinx重合。
首先用mesh函数.
2;
k=[0.2,0.1,0.05];
[xi,yi]=meshgrid(-10:
k(1,i+1):
10);
z=exp(-abs(xi))+cos(xi+yi)+1./(xi.^2+yi.^2+1);
subplot(2,2,i+1)
title('
mesh'
mesh(xi,yi,z)
再用surf函数.
第二章实验
1.列主元三角分解A.
clc
A=[11111;
12345;
1361015;
14102035;
15153570]
E=eye(5)
[m,n]=size(A);
ifm~=n
error('
不是方阵'
)
return
end
ifdet(A)==0
不能分解'
u=A;
p=eye(m);
l=eye(m);
fori=1:
m
forj=i:
t(j)=u(j,i);
fork=1:
i-1
t(j)=t(j)-u(j,k)*u(k,i);
a=i;
b=abs(t(i));
forj=i+1:
ifb<
abs(t(j))
b=abs(t(j));
a=j;
ifa~=i
forj=1:
m
c=u(i,j);
u(i,j)=u(a,j);
u(a,j)=c;
c=p(i,j);
p(i,j)=p(a,j);
p(a,j)=c;
c=t(a);
t(a)=t(i);
t(i)=c;
u(i,i)=t(i);
u(j,i)=t(j)/t(i);
i-1
u(i,j)=u(i,j)-u(i,k)*u(k,j);
l=tril(u,-1)+eye(m)
u=triu(u,0)
B=A/E
2.追赶法求方程,电路的电流
A=[2-2000000;
-25-200000;
0-25-20000;
00-25-2000;
000-25200;
0000-25-20;
00000-25-2;
000000-25]
bb=[220/270000000]
n=size(A,1);
s=zeros(n,1);
%-----È
¡
³
ö
È
ý
¶
Ô
½
Ç
--------
b=diag(A);
a=diag(A,-1);
c=diag(A,1);
d=zeros(n,1);
u=zeros(n-1,1);
fori=1:
n-1
d
(1)=b
(1);
u(i)=c(i)/d(i);
d(i+1)=b(i+1)-a(i)*u(i);
%-----×
·
µ
Ä
¹
Ì
------------
y=zeros(n,1);
y
(1)=bb
(1)/d
(1);
fori=2:
n
y(i)=(bb(i)-a(i-1)*y(i-1))/d(i);
%-----¸
Ï
---------------
s(n)=y(n);
fori=n-1:
-1:
1
s(i)=y(i)-u(i)*s(i+1);
s
3.方程组的性态与矩阵条件数的实验
clc
n=5;
A=zeros(n);
C=zeros(n);
b=zeros(1,n);
x(i)=1+0.1*i;
A(i,j)=x(i)^(j-1);
C(i,j)=1/(i+j-1);
b(i)=b(i)+A(i,j);
A
C
b
d=cond(A,2)
DD=A\b'
(1)当n=5,10,20时,cond(A)=5.3615e+05,8.6823e+11,1.5428e+22。
看出越来越病态了。
(2)当n=5,10,20时,解分别如下:
说明条件数越大,越病态,发散了。
当n=10时就一个值有一定误差。
(3)n=10,解方程(A+delta)x=b.
clc
n=10;
%C(i,j)=1/(i+j-1);
A(2,2)=A(2,2)+10^(-8)
A(10,10)=A(10,10)+10^(-8)
x=A\b'
那么结果如下:
可以看出很小的扰动导致解误差变得较大了。
误差
(1)行列式,条件数,特征值
(2)(A+A0)(x+x0)=b,求此时的x0与||x0||2.
代码如下:
A=[10787;
7565;
86109;
75910]
A1=[107.28.16.9;
7.085.076.025;
8.25.899.969.01;
6.985.048.979.98]
b=[32233331]
A0=A1-A
x1=A1\b'
x0=x1-x
fanshux0=sqrt(x0(1,1)^2+x0(2,1)^2+x0(3,1)^2+x0(4,1)^2)
fanshux=sqrt(1+1+1+1)
eigAA=eig(A'
*A);
fanshuA=sqrt(max(eigAA))
eigA0A0=eig(A0'
*A0);
fanshuA0=sqrt(max(eigA0A0))
x0_x=fanshux0/fanshux
A0_A=fanshuA0/fanshuA
condA=cond(A)
x0_xx=(condA*A0_A)/(1-condA*A0_A)
结果如下:
然后比较:
然后利用公式计算:
左边=0.8225,右边=1.0358
则有0.8225<
1.0358。
第三章实验
1.分别用Jacobi,Gauss-Seild,共轭梯度法解方程
(1)Jacobi迭代法
A=[101234;
19-12-3;
2-173-5;
32312-1;
4-3-5-115]
b=[12-2714-1712]
n=length(A);
k=0;
ep=10^(-7)
it_max=100;
x=zeros(n,1);
y=zeros(n,1);
fprintf('
x1x2x3x4x5\n'
while1
n
y(i)=b(i);
ifj~=i
y(i)=y(i)-A(i,j)*x(j);
ifabs(A(i,i))<
1e-10|k==it_max
return;
y(i)=y(i)/A(i,i);
ifnorm(y-x,inf)<
ep
x=y;
k=k+1;
%d,%f%f%f%f%f\n'
k,x(1,1),x(2,1),x(3,1),x(4,1),x(5,1))
得的结果为
(2)Gauss-Seild迭代法
x0=[00000]'
;
ep=10^(-7);
N=100;
k=1;
\n**********start**************\n'
k);
fork=1:
N
%d\t'
%calculatex
(1)
x
(1)=b
(1)/A(1,1);
forj=2:
x
(1)=x
(1)-A(1,j)*x0(j)/A(1,1);
%f\t'
x
(1));
%calculatex(i)
fori=2:
x(i)=b(i)/A(i,i);
x(i)=x(i)-A(i,j)*x(j)/A(i,i);
x(i)=x(i)-A(i,j)*x0(j)/A(i,i);
x(i));
%calculatex(n)
x(n)=b(n)/A(n,n);
x(n)=x(n)-A(n,j)*x(j)/A(n,n);
x(n));
ifnorm(x-x0,inf)<
else
x0=x;
\n'
);
\n***********end*************\n'
x0
运行结果如下:
可以看出Gauss-Seild只要41次就可以计算结果,而Jacobi需要76次。
(3)共轭梯度迭代法
b=[12-2714-1712]'
N=length(A);
eps=10^(-7);
x=zeros(N,1);
%µ
ü
´
ú
Ë
Æ
ò
Á
¿
.normÊ
¾
Ø
Õ
ó
£
¬
Ê
r=b-A*x;
%Æ
ä
±
í
¸
º
Ñ
÷
fai=1/2x'
AX-bx
d=r;
fork=0:
N-1
Ú
%d´
Î
'
k+1);
a=(norm(r)^2)/(d'
*A*d)
x=x+a*d
rr=b-A*x;
%rr=r(k+1)
if(norm(rr)<
=eps)||(k==N-1)
B=(norm(rr)^2)/(norm(r)^2);
d=rr+B*d;
r=rr;
计算结果如下:
理论上只需要N次得到精确值
2.利用共轭梯度法计算矩阵--10^5阶
n=10^5;
ifi==j
A(i,j)=3;
elseifabs(i-j)==1
A(i,j)=-1;
elseif(i+j==n+1)&
&
(i~=n/2)&
(i~=n/2+1)
A(i,j)=1/2;
A(i,j)=0;
b(1,1)=2.5;
b(1,2:
10^5/2-1)=1.5;
b(1,10^5/2)=1.0;
b(1,10^5/2+1)=1.0;
b(1,10^5/2+2:
10^5-1)=1.5;
b(1,10^5)=2.5;
b=b'
fai=1/2*x'
A=(norm(r)^2)/(d'
*A*d);
x=x+A*d;
x'
3.利用cgs,bicg,bicgstab,等计算矩阵解
利用cgs()函数
A=gallery('
wilk'
21)
b=sum(A,2);
tol=1e-12;
maxit=15;
M1=diag([10:
111:
10]);
x=cgs(A,b,tol,maxit,M1)
利用bicg()函数
x=bicg(A,b,tol,maxit,M1)
利用bicgstab()函数
maxit=10;
x=bicgstab(A,b,tol,maxit,M1)
第五章实验
1.Newton插值f(x)=1/(1+4x^2).图形,误差
(1)输出插值多项式。
程序:
function[p]=NewtonChazhi(x,y,n)
%UNTITLED3此处显示有关此函数的摘要
%此处显示详细说明
f=zeros(n,n);
f(:
1)=y'
forj=2:
fori=j:
f(i,j)=(f(i,j-1)-f(i-1,j-1))/(x(i)-x(i-j+1));
p=f(n,n);
fork=n:
p=conv(p,poly(x(k-1)));
d=length(p);
p(d)=p(d)+f(k-1,k-1);
p
disp('
牛顿插值多项式:
polynomial=poly2sym(p,'
x'
x=-5:
5;
y=1./(1+4*x.^2);
n=length(x);
[p]=NewtonChazhi(x,y,n)
结果:
polynomial=
-(7319042784910035*x^10)/147573952589676412928+(256*x^8)/93425+(5*x^7)/1152921504606846976-(7410620163505401*x^6)/144115188075855872+x^5/9007199254740992+(36624*x^4)/93425+x^3/9007199254740992-(5148893614132309*x^2)/4503599627370496+(19*x)/36028797018963968+1
p=
Columns1through6
-0.000000.00270.0000-0.05140.0000
Columns7through11
1.39200.0000-1.14330.00001.0000
(2)在[-5,5]均匀插入99个节点,计算这些节点上函数f(x)的近似值,在同一张图上画出原函数与插值函数的图形。
x=-5:
X=linspace(-5,5,101);
y1=1./(1+4*X.^2);
101
N(i)=polyval(p,X(i));
E(i)=N(i)-y1(i);
plot(X,y1,X,N)
(4)误差图形
2.f(x)=1/(1+4x^2)插值
3.飞机的外形轮廓
y=[0,-30,-36,-35,-28.44,-9.4,0];
x2=520:
-8.67:
0;
x=[520,280,156.6,78,39.62,3.1,0];
y2=spline(x,y,x2);
%用样条插值计算近似值
x1=[0,3.1,39.62,78,156.6,280,520];
y1=[0,9.4,28.44,35,36,30,0];
x3=0:
8.67:
520;
y3=spline(x1,y1,x3);
plot(x2,y2,'
r-'
x3,y3,'
g-'
title('
机翼外形曲线'
第九章四阶龙格库塔方法(自选)
四阶龙格-库塔(Runge-Kutta)方法,因为在本人研究方向上,要计算姿态结算方程组,而经典的方法就是四阶RK方法,故这章自己选题做这个方面的实验。
1.解算微分方程组
导航姿态方程如下:
角度与角速度关系
2.Matlab的四阶Rk方法
3.作图比较,一定误差
黑线表示实验测试计算结果,红线绿线表示matlab计算结果。
两者误差较大。
由于计算中全是三角函数计算,误差大。