数值分析课程设计.docx

数值分析课程设计.docx

《数值分析课程设计.docx》由会员分享，可在线阅读，更多相关《数值分析课程设计.docx（17页珍藏版）》请在冰点文库上搜索。

数值分析课程设计

1

2用高斯消元法的消元过程作矩阵分解。

设

消元过程可将矩阵A化为上三角矩阵U，试求出消元过程所用的乘数

、

、

并以如下格式构造下三角矩阵L和上三角矩阵U

验证：

矩阵A可以分解为L和U的乘积，即A=LU。

Solution:

function[L,U,flag]=LU_Decom（A）

[n,m]=size（A）;

ifn~=m

error（'TherowsandcolumnsofmatrxAmustbeequal!

'）;

return;

end

L=eye（n）;U=zeros（n）;flag='OK';

fork=1:

n

forj=k:

n

z=0;

forq=1:

k-1

z=z+L（k,q）*U（q,j）;

end

U（k,j）=A（k,j）-z;

end

ifabs（U（k,k））

flag='failure';return;

end

fori=k+1:

n

z=0;

forq=1:

k-1

z=z+L（i,q）*U（k,k）;

end

L（i,q）=（A（i,k）-z）/U（k,k）;

end

end

输入命令：

>>A=[20,2,3;1,8,1;2,-3,15];

>>[L,U,flag]=LU_Decom（A）

运行结果：

L=

1.000000

01.00000

-0.375001.0000

U=

20.00002.00003.0000

08.00001.0000

0016.1250

flag=

OK

>>

3

4

5

6．用欧拉公式和四阶龙格-库塔法分别求解下列初值问题；

（1）

欧拉方法

functionE=Euler_1（fun,x0,y0,xN,N）

x=zeros（1,N+1）;y=zeros（1,N+1）;

x

（1）=x0;y

（1）=y0;

h=（xN-x0）/N;

forn=1:

N

x（n+1）=x（n）+h;

y（n+1）=y（n）+h*feval（fun,x（n）,y（n））;

end

T=[x',y']

functionz=f（x,y）

z=0.9*y/（1+2*x）;

输入命令

>>Euler_1（'f',0,1,1,10）

运行结果

T=

01.0000

0.10001.0900

0.20001.1718

0.30001.2471

0.40001.3172

0.50001.3831

0.60001.4453

0.70001.5045

0.80001.5609

0.90001.6149

1.00001.6668

>>

四阶龙格-库塔方法

functionR=rk4（f,a,b,ya,N）

h=（b-a）/N;

T=zeros（1,N+1）;

Y=zeros（1,N+1）;

T=a:

h:

b;

Y

（1）=ya;

forj=1:

N

k1=h*feval（f,T（j）,Y（j））;

k2=h*feval（f,T（j）+h/2,Y（j）+k1/2）;

k3=h*feval（f,T（j）+h/2,Y（j）+k2/2）;

k4=h*feval（f,T（j）+h,Y（j）+k3）;

Y（j+1）=Y（j）+（k1+2*k2+2*k3+k4）/6;

end

R=[T'Y']

functionz=（x,y）

z=0.9*y/（1+2*x）

输入命令

>>rk4（'f',0,1,1,10）

运行结果

R=

01.0000

0.10001.0855

0.20001.1635

0.30001.2355

0.40001.3028

0.50001.3660

0.60001.4259

0.70001.4828

0.80001.5372

0.90001.5894

1.00001.6395

ans=

01.0000

0.10001.0855

0.20001.1635

0.30001.2355

0.40001.3028

0.50001.3660

0.60001.4259

0.70001.4828

0.80001.5372

0.90001.5894

1.00001.6395

>>

（2）欧拉方法

functionE=Euler_1（fun,x0,y0,xN,N）

x=zeros（1,N+1）;y=zeros（1,N+1）;

x

（1）=x0;y

（1）=y0;

h=（xN-x0）/N;

forn=1:

N

x（n+1）=x（n）+h;

y（n+1）=y（n）+h*feval（fun,x（n）,y（n））;

end

T=[x',y']

functionz=（x,y）

z=-x*y/（1+x^2）;

输入命令

>>Euler_1（'f',0,1,1,10）

运行结果

T=

01.0000

0.10001.0900

0.20001.1718

0.30001.2471

0.40001.3172

0.50001.3831

0.60001.4453

0.70001.5045

0.80001.5609

0.90001.6149

1.00001.6668

>>

四阶龙格-库塔方法

functionR=rk4（f,a,b,ya,N）

h=（b-a）/N;

T=zeros（1,N+1）;

Y=zeros（1,N+1）;

T=a:

h:

b;

Y

（1）=ya;

forj=1:

N

k1=h*feval（f,T（j）,Y（j））;

k2=h*feval（f,T（j）+h/2,Y（j）+k1/2）;

k3=h*feval（f,T（j）+h/2,Y（j）+k2/2）;

k4=h*feval（f,T（j）+h,Y（j）+k3）;

Y（j+1）=Y（j）+（k1+2*k2+2*k3+k4）/6;

end

R=[T'Y']

functionz=（x,y）

z=-x*y/（1+x^2）;

输入命令

>>rk4（'f',0,1,1,10）

运行结果

R=

01.0000

0.10001.0855

0.20001.1635

0.30001.2355

0.40001.3028

0.50001.3660

0.60001.4259

0.70001.4828

0.80001.5372

0.90001.5894

1.00001.6395

ans=

01.0000

0.10001.0855

0.20001.1635

0.30001.2355

0.40001.3028

0.50001.3660

0.60001.4259

0.70001.4828

0.80001.5372

0.90001.5894

1.00001.6395

>>

7用二分法求下列方程在指定区间[a，b]上的实根近似值：

（要求误差不超过0.01）

（1）x-sinx-1=0，[a,b]=[1,3]；

（2）xsinx=1，[a,b]=[0,1.5].

（1）functionroot=bisect（f,a,b,eps）

if（nargin==3）

eps=1.0e-4;

end

f1=subs（sym（f）,findsym（sym（f））,a）;

f2=subs（sym（f）,findsym（sym（f））,b）;

if（f1==0）

root=a;

end

if（f2==0）

root=b;

end

if（f1*f2>0）

disp（'Á½¶Ëµãº¯ÊýÖµ³Ë»ý´óÓÚ0!

'）;

return;

else

root=FindRoots（f,a,b,eps）;

end

functionr=FindRoots（f,a,b,eps）

f_1=subs（sym（f）,findsym（sym（f））,a）;

f_2=subs（sym（f）,findsym（sym（f））,b）;

mf=subs（sym（f）,findsym（sym（f））,（a+b）/2）;

if（f_1*mf>0）

t=（a+b）/2;

r=FindRoots（f,t,b,eps）;

else

if（f_1*mf==0）

r=（a+b）/2;

else

if（abs（b-a）<=eps）

r=（b+3*a）/4;

else

s=（a+b）/2;

r=FindRoots（f,a,s,eps）;

end

end

end

输入命令

>>root=bisect（'x-sin（x）-1',1,3）

运行结果

root=

1.9346

（2）functionroot=bisect（f,a,b,eps）

if（nargin==3）

eps=1.0e-4;

end

f1=subs（sym（f）,findsym（sym（f））,a）;

f2=subs（sym（f）,findsym（sym（f））,b）;

if（f1==0）

root=a;

end

if（f2==0）

root=b;

end

if（f1*f2>0）

disp（'Á½¶Ëµãº¯ÊýÖµ³Ë»ý´óÓÚ0!

'）;

return;

else

root=FindRoots（f,a,b,eps）;

end

functionr=FindRoots（f,a,b,eps）

f_1=subs（sym（f）,findsym（sym（f））,a）;

f_2=subs（sym（f）,findsym（sym（f））,b）;

mf=subs（sym（f）,findsym（sym（f））,（a+b）/2）;

if（f_1*mf>0）

t=（a+b）/2;

r=FindRoots（f,t,b,eps）;

else

if（f_1*mf==0）

r=（a+b）/2;

else

if（abs（b-a）<=eps）

r=（b+3*a）/4;

else

s=（a+b）/2;

r=FindRoots（f,a,s,eps）;

end

end

end

输入命令

>>root=bisect（'x*sin（x）-1',0,1.5）

运行结果

root=

1.1142

8分别用直接法、雅可比迭代法、赛德尔迭代法求解线性方程组AX=b。

（1）

，

（2）

，

（1）function[x,n]=jacobi（A,b,x0,eps,varargin）

ifnargin==3

eps=1.0e-6;

M=200;

elseifnargin<3

error

return

elseifnargin==5

M=varargin{1};

end

D=diag（diag（A））;

L=-tril（A,-1）;

U=-triu（A,1）;

B=D\（L+U）;

f=D\b;

x=B*x0+f;

n=1;

whilenorm（x-x0）>=eps

x0=x;

x=B*x0+f;

n=n+1;

if（n>=M）

disp（'Warning:

µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²!

'）;

end

end

输入系数矩阵和右端项

A=[4100000000;1410000000;0000000000;0000000000;0000000000;0000000000;0000000000;0000000000;0000000141;0000000014],

A=

4100000000

1410000000

0000000000

0000000000

0000000000

0000000000

0000000000

0000000000

0000000141

0000000014

>>b=[2100000012]'

b=

2

1

0

0

0

0

0

0

1

2

输入命令

>>x0=ones（10,1）;

>>[x,n]=Jacobi（A,b,x0）

输出结果

x=

0.2500

-0.2500

NaN

NaN

NaN

NaN

NaN

NaN

NaN

NaN

n=

1

（2）function[x,n]=jacobi（A,b,x0,eps,varargin）

ifnargin==3

eps=1.0e-6;

M=200;

elseifnargin<3

error

return

elseifnargin==5

M=varargin{1};

end

D=diag（diag（A））;

L=-tril（A,-1）;

U=-triu（A,1）;

B=D\（L+U）;

f=D\b;

x=B*x0+f;

n=1;

whilenorm（x-x0）>=eps

x0=x;

x=B*x0+f;

n=n+1;

if（n>=M）

disp（'Warning:

µü´ú´ÎÊýÌ«¶à£¬¿ÉÄܲ»ÊÕÁ²!

'）;

end

end

输入系数矩阵和右端项

9,

10,