数值代数上机实验报告.docx

资源描述

数值代数上机实验报告.docx

《数值代数上机实验报告.docx》由会员分享，可在线阅读，更多相关《数值代数上机实验报告.docx（22页珍藏版）》请在冰豆网上搜索。

数值代数上机实验报告.docx

数值代数上机实验报告

数值代数课程设计实验报告

姓名：

班级：

学号：

实验日期：

一、实验名称

代数的数值解法

二、实验环境

MATLAB7.0

实验一、平方根法与改进平方根法

一、实验要求：

用熟悉的计算机语言将不选主元和列主元Gasuss消元法编写成通用的子程序，然后用编写的程序求解下列方程组

用所编的程序分别求解40、84、120阶方程组的解。

2、算法描述及实验步骤

GAuss程序如下：

（1）求A的三角分解：

;

（2）求解

得y；

（3）求解

得x；

列主元Gasuss消元法程序如下：

1求A的列主元分解：

;

2求解

得y；

3求解

得x；

三、调试过程及实验结果：

%----------------方程系数----------------

>>A1=Sanduijiaozhen（8,6,1,40）;

>>A2=Sanduijiaozhen（8,6,1,84）;

>>A3=Sanduijiaozhen（8,6,1,120）;

>>b1

（1）=7;b2

（1）=7;b3

（1）=7;

>>fori=2:

b1（i）=15;

end

>>b1（40）=14;

>>fori=2:

b2（i）=15;

end

>>b2（40）=14;

>>fori=2:

119

b1（i）=15;

end

>>b3（120）=14;

%----------------方程解----------------

>>x11=GAuss（A1,b1'）

>>x12=GAussZhu（A1,b1'）

>>x21=GAuss（A2,b2'）

>>x22=GAussZhu（A3,b3'）

>>x31=GAuss（A3,b3'）

>>x32=GAussZhu（A3,b3'）

运行结果：

（n=40）

GAuss消元法的解即为

x11=

1.0000

列主元GAuss消元法的解即为

x12=

6、源程序：

functionA=Sanduijiaozhen（a,b,c,n）

%生成n阶以a,b,c为元素的三对角阵

A=diag（b*ones（1,n）,0）+diag（c*ones（1,n-1）,1）+diag（a*ones（1,n-1）,-1）;

functionx=GAuss（A,b）

n=length（b）;

x=b;

%-------分解---------------

fori=1:

n-1

forj=i+1:

mi=A（j,i）/A（i,i）;

b（j）=b（j）-mi*b（i）;

fork=i:

A（j,k）=A（j,k）-mi*A（i,k）;

end

AB=[A,b];

end

%-----------回代------------------

x（n）=b（n）/A（n,n）;

fori=n-1:

-1:

s=0;

forj=i+1:

s=s+A（i,j）*x（j）;

end

x（i）=（b（i）-s）/A（i,i）;

end

functionx=GAussZhu（A,b）

n=length（b）;

x=b;

%----------------------选主元---------------------

fork=1:

n-1

a_max=0;

fori=k:

ifabs（A（i,k））>a_max

a_max=abs（A（i,k））;r=i;

end

ifr>k

forj=k:

z=A（k,j）;A（k,j）=A（r,j）;A（r,j）=z;

end

z=b（k）;b（k）=b（r）;b（r）=z;

end

%--------------消元-----------------

fori=k+1:

m=A（i,k）/A（k,k）;

forj=k:

A（i,j）=A（i,j）-m*A（k,j）;

end

b（i）=b（i）-m*b（k）;

end

ifabs（A（n,n））==0

return;

end

AbZhu=[A,b];

%----------------回代-----------------------

x（n）=b（n）/A（n,n）;

fori=n-1:

-1:

forj=i+1:

b（i）=b（i）-A（i,j）*x（j）;

end

x（i）=b（i）/A（i,i）;

end

实验二、平方根法与改进平方根法

一、实验要求：

用计算机语言将平方根法和改进的平方根法编成通用的子程序，然后用编写的程序求解对称正定方程组

阶方程组AX=b，

二、算法描述及实验步骤：

平方根法函数程序如下：

1、求A的Cholesky分解：

;

2、求解

得y；

3、求解

得x；

改进平方根法函数程序如下：

1、求A的Cholesky分解：

;

2、求解

得y；

3、求解

得x；

三、调试过程及实验结果：

clear;clc;

%----------------方程系数----------------

>>A=Sanduijiaozhen（1,10,1,100）;

>>b

（1）=11;

>>fori=2:

b（i）=12;

end

>>b（100）=11;

>>x1=Cholesky（A,b'）

>>x2=GJCholesky（A,b'）

运行结果：

平方根法的解即为

x1=

1.0000

改进平方根法解得的解即为

x2=

1.0000

0.9999

1.0009

0.9908

1.0908

0.1010

四、源程序：

functionx=Cholesky（A,b）

n=size（A）;

n=n

（1）;

%x=A^-1*b;

%disp（'Matlab自带解即为x'）;

%-----------------Cholesky分解-------------------

fork=1:

A（k,k）=sqrt（A（k,k））;

A（k+1:

n,k）=A（k+1:

n,k）/A（k,k）;

forj=k+1:

A（j:

n,j）=A（j:

n,j）-A（j:

n,k）*A（j,k）;

end

%------------------前代法求解Ly=b----------------------------

forj=1:

n-1

b（j）=b（j）/A（j,j）;

b（j+1:

n）=b（j+1:

n）-b（j）*A（j+1:

n,j）;

end

b（n）=b（n）/A（n,n）;

%-----------------回代法求解L'x=y-----------------------------

A=A';

forj=n:

-1:

b（j）=b（j）/A（j,j）;

b（1:

j-1）=b（1:

j-1）-b（j）*A（1:

j-1,j）;

end

（1）=b

（1）/A（1,1）;

disp（'平方根法的解即为'）;

functionb=GJCholesky（A,b）

n=size（A）;

n=n

（1）;

v=zeros（n,1）;

%----------------------LDL'分解-----------------------------

forj=1:

fori=1:

j-1

v（i）=A（j,i）*A（i,i）;

end

A（j,j）=A（j,j）-A（j,1:

j-1）*v（1:

j-1）;

A（j+1:

n,j）=（A（j+1:

n,j）-A（j+1:

n,1:

j-1）*v（1:

j-1））/A（j,j）;

end

B=diag（A）;

D=zeros（n）;

fori=1:

D（i,i）=B（i）;

A（i,i）=1;

end

%-------------------前代法---------------------------

A=tril（A）;%得到L和D

forj=1:

n-1

b（j）=b（j）/A（j,j）;

b（j+1:

n）=b（j+1:

n）-b（j）*A（j+1:

n,j）;

end

b（n）=b（n）/A（n,n）;

%-----------------回代法-----------------------------

A=D*（A'）;

forj=n:

-1:

b（j）=b（j）/A（j,j）;

b（1:

j-1）=b（1:

j-1）-b（j）*A（1:

j-1,j）;

end

（1）=b

（1）/A（1,1）;

disp（'改进平方根法解得的解即为'）;

实验三、二次多项式拟合

一、实验要求：

用计算机语言编制利用QR分解求解线性最小二乘问题的通用子程序，用编写的程序求解一个二次多项式

使在残向量的

范数最小的意义下拟合下面的数据

-1-0.75-0.500.250.50.75

1.000.81250.751.001.31251.752.3125

二、算法描述及实验步骤：

QR分解求解程序如下：

1、求A的

分解;

2、计算

;

3、求解上三角方程

得x；

三、调试过程及实验结果：

>>t=[-1-0.75-0.500.250.50.75];

>>y=[1.000.81250.751.001.31251.752.3125];

>>plot（t,y,'r*'）;

>>legend（'实验数据（ti,yi）'）;

>>xlabel（'t'）,ylabel（'y'）;

>>title（'二次多项式拟合的数据点（ti,yi）的散点图'）;

运行后屏幕显示数据的散点图（略）.

（3）编写下列MATLAB程序计算

在

处的函数值，即输入程序

>>symsabc

>>t=[-1-0.75-0.500.250.50.75];

>>fi=a.*t.^2+b.*t+c

%运行后屏幕显示关于

的线性方程组

fi=

[a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b+c]

编写构造残向量2范数的MATLAB程序

>>y=[1.000.81250.751.001.31251.752.3125];

>>fy=fi-y;fy2=fy.^2;J=sum（fy.^2）;

运行后屏幕显示误差平方和如下

（a-b+c-1）^2+（9/16*a-3/4*b+c-13/16）^2+（1/4*a-1/2*b+c-3/4）^2+（c-1）^2+（1/16*a+1/4*b+c-21/16）^2+（1/4*a+1/2*b+c-7/4）^2+（9/16*a+3/4*b+c-37/16）^2

为求

使

达到最小，只需利用极值的必要条件

，

，得到关于

的线性方程组，这可以由下面的MATLAB程序完成，即输入程序

>>Ja1=diff（J,a）;Ja2=diff（J,b）;Ja3=diff（J,c）;

>>Ja11=simple（Ja1）,Ja21=simple（Ja2）,Ja31=simple（Ja3）

运行后屏幕显示J分别对

的偏导数如下

Ja11=

451/128*a-63/32*b+43/8*c-887/128

Ja21=

-63/32*a+43/8*b-3/2*c-61/32

Ja31=

43/8*a-3/2*b+14*c-143/8

解线性方程组

，输入下列程序

>>A=[451/128,-63/32,-3/2;-63/32,43/8,-3/2;43/8,-3/2,14];

>>B=[887/128,61/32,143/8];

>>C=B/A,f=poly2sym（C）

运行后屏幕显示拟合函数f及其系数C如下

0.30810.85871.4018

924/2999*x^2+10301/11996*x+4204/2999

故所求的拟合曲线为

四、源程序：

>>t=[-1-0.75-0.500.250.50.75];

>>y=[1.000.81250.751.001.31251.752.3125];

>>plot（t,y,'r*'）;

>>legend（'实验数据（ti,yi）'）;

>>xlabel（'t'）,ylabel（'y'）;

>>title（'二次多项式拟合的数据点（ti,yi）的散点图'）;

>>symsabc

>>t=[-1-0.75-0.500.250.50.75];

>>fi=a.*t.^2+b.*t+c

fi=

[a-b+c,9/16*a-3/4*b+c,1/4*a-1/2*b+c,c,1/16*a+1/4*b+c,1/4*a+1/2*b+c,9/16*a+3/4*b+c]

>>y=[1.000.81250.751.001.31251.752.3125];

>>fy=fi-y;fy2=fy.^2;J=sum（fy.^2）

（a-b+c-1）^2+（9/16*a-3/4*b+c-13/16）^2+（1/4*a-1/2*b+c-3/4）^2+（c-1）^2+（1/16*a+1/4*b+c-21/16）^2+（1/4*a+1/2*b+c-7/4）^2+（9/16*a+3/4*b+c-37/16）^2

>>Ja1=diff（J,a）;Ja2=diff（J,b）;Ja3=diff（J,c）;

>>Ja11=simple（Ja1）,Ja21=simple（Ja2）,Ja31=simple（Ja3）

Ja11=

451/128*a-63/32*b+43/8*c-887/128

Ja21=

-63/32*a+43/8*b-3/2*c-61/32

Ja31=

43/8*a-3/2*b+14*c-143/8

>>A=[451/128,-63/32,-3/2;-63/32,43/8,-3/2;43/8,-3/2,14];

>>B=[887/128,61/32,143/8];

>>C=B/A,f=poly2sym（C）

0.30810.85871.4018

924/2999*x^2+10301/11996*x+4204/2999

展开阅读全文