数值分析大作业

格式：docx
大小：71.90 KB
文档页数：14

数值分析大作业

姓名：黄晨晨

学号：S1*******

学院：储运与建筑工程学院学院班级：储建研17-2

实验3.1 Gauss消去法的数值稳定性实验

实验目的：

理解高斯消元过程中出现小主元即很小时引起方程组解数值不定性实验内容：

求解方程组Ax=b，其中

（1）A1=0.3×10?1559.1431

5.291?

6.130?12

11.2952

1211

，b1=

59.17

46.78

;

（2）A2=10?701

3 2.09999999999962

5?15?1

0102

，b2=

5.900000000001

5 1

;

实验要求：

(1)计算矩阵的条件数，判断系数矩阵是良态的还是病态的

(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4

(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4

(4)观察小主元并分析对计算结果的影响

（1）计算矩阵的条件数，判断系数矩阵是良态的还是病态的代码：

format longe

format compact

A1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1]

b1=[59.17;46.78;1;2]

n=4

C1=cond(A1,1) %C1为A1矩阵1范数下的条件数

C2=cond(A1,2) %C2为A1矩阵2范数下的条件数

C3=cond(A1,inf) %C3为1矩阵谱范数下的条件数

结果：

C1 =

1.362944708720448e+02

C2 =

6.842955771253409e+01

C3 =

8.431146*********e+01

显然A1矩阵为病态矩阵

将矩阵A2，b2输入上述代码中求得A2矩阵的条件数为：C1 =

1.928316831682894e+01

C2 =

8.993938090170119e+00

C3 =

1.835643564356072e+01 显然A2矩阵也为病态矩阵

(2)用Gauss列主元消去法求得L和U及解向量x1,x2∈R4

代码：

for k=1:n-1

a=max(abs(A1(k:n,k)))

if a==0

return

end

for i=k:n

if abs(A1(i,k))==a

y=A1(i,:)

A1(i,:)=A1(k,:)

A1(k,:)=y

x=b1(i,:)

b1(i,:)=b1(k,:)

b1(k,:)=x

break

end

if A1(k,k)~=0

A1(k+1:n,k)=A1(k+1:n,k)/A1(k,k)

A1(k+1:n,k+1:n)=A1(k+1:n,k+1:n)-A1(k+1:n,k)*A1(k,k+1:n)

else

break

end

L=tril(A1,0);

for i=1:n

L(i,i)=1;

end L

U=triu(A1,0)

y1=L\b1

x1=U\y1

得到如下结果：

x1 =3.845714853511634e+00

1.609517394778522e+00

-1.547605454206655e+01

1.041130489899787e+01

将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2]

b2=[8;5.900000000001;5;1]

代入上述代码求得结果如下：

X2 =

4.440892098500626e-16

-9.999999999999993e-01

9.999999999999997e-01

1.000000000000000e+00

(3)用不选主元的高斯消去法求得L和U及解向量x1,x2∈R4代码： format longe format compact

A1=[0.3*10^-15,59.14,3,1;5.291,-6.130,-1,2;11.2,9,5,2;1,2,1,1]

b1=[59.17;46.78;1;2]

[L,U]=lu(A1)

y1=L\b1

x1=U\y1

求得如下结果：

x1=

3.845714853511634e+00

1.609517394778522e+00

-1.547605454206655e+01

1.041130489899787e+01

将A2=[10,-7,0,1;-3,2.0999********,6,2;5,-1,5,-1;0,1,0,2]

b2=[8;5.900000000001;5;1]

代入上述代码，求得结果如下：

x 2 =

4.440892098500626e-16 -9.999999999999993e-01

9.999999999999997e-01 9.999999999999999e-01

（2）（3）求得结果相同，可验证结果正确。 (4)观察小主元并分析对计算结果的影响实验结果表明，小主元的存在对计算结果影响较大。小主元的存在会使得计算结果有较大的误差。实验3.2方程组的性态和条件数实验实验目的：

理解条件数的意义和方程组的性态对解向量的影响。实验内容：

已知两个方程组A1x=b 和A2x=b ，其中：

A1= 1

x 0x 02…x 0n

1x 1x 12...x 1n 1x 2x 22...x 2n (1)

xn 2…xn n

A2= 11/21/3…1/n

1/21/31/4…1/(n +1)

1/3

1/41/5…1/(n +2)……………1/n 1/(n +1)1/(n +2)

…

1/(2n ?1)

b= a 1j n j =1, a 2j n j =1,… a 1nj n j =1, T

实验要求：

对A1，取xk=1+0.1k,k=0,1,…，n，下面均用Matlab函数“x=A/b”计算方程组的解。

（1）取n=4,6,8,分别求出A1,A2的条件数，判别它们是否是病态阵？随n的增大，矩阵性态的变化如何？

代码：

format shorte

format compact

A1=[1,1,1,1,1;1,1.1,1.1^2,1.1^3,1.1^4;1,1.2,1.2^2,1.2^3,1.2^4;1,1.3,1.3^2 ,1.3^3,1.3^4;1,1.4,1.4^2,1.4^3,1.4^4]

C1=cond(A1,1) %C1为A1矩阵1范数下的条件数

C2=cond(A1,2) %C2为A1矩阵2范数下的条件数 C3=cond(A1,inf) %C3为1矩阵谱范数下的条件数

1.对A1矩阵运算

n=4时：

C1 =

6.5120e+05

C2 =

3.5740e+05

C3 =

6.2755e+05

n=4时，A1为病态矩阵

n=6时：

A1=[1,1,1,1,1,1,1;1,1.1,1.1^2,1.1^3,1.1^4,1.1^5,1.1^6;1,1.2,1.2^2,1.2^3,1 .2^4,1.2^5,1.2^6;1,1.3,1.3^2,1.3^3,1.3^4,1.3^5,1.3^6;1,1.4,1.4^2,1.4^3,1.

4^4,1.4^5,1.4^6;1,1.5,1.5^2,1.5^3,1.5^4,1.5^5,1.5^6;1,1.6,1.6^2,1.6^3,1. 6^4,1.6^5,1.6^6]

C1 =

1.8531e+08

C2 =

8.7385e+07

C3 =

1.6574e+08

n=6时，A1为病态矩阵

n=8时：

A1=[1,1,1,1,1,1,1,1,1;1,1.1,1.1^2,1.1^3,1.1^4,1.1^5,1.1^6,1.1^7,1.1^8;1,1 .2,1.2^2,1.2^3,1.2^4,1.2^5,1.2^6,1.2^7,1.2^8;1,1.3,1.3^2,1.3^3,1.3^4,1.3

^5,1.3^6,1.3^7,1.3^8;1,1.4,1.4^2,1.4^3,1.4^4,1.4^5,1.4^6,1.4^7,1.4^8;1,1 .5,1.5^2,1.5^3,1.5^4,1.5^5,1.5^6,1.5^7,1.5^8;1,1.6,1.6^2,1.6^3,1.6^4,1.6 ^5,1.6^6,1.6^7,1.6^8;1,1.7,1.7^2,1.7^3,1.7^4,1.7^5,1.7^6,1.7^7,1.7^8;1,1 .8,1.8^2,1.8^3,1.8^4,1.8^5,1.8^6,1.8^7,1.8^8]

C1 =

5.0565e+10

C2 =

2.2739e+10

C3 =

4.4788e+10

n=8时，A1为病态矩阵

2.对A2矩阵运算

n=4时：

A2=[1,1/2,1/3,1/4;1/2,1/3,1/4,1/5;1/3,1/4,1/5,1/6;1/4,1/5,1/6,1/7]

C1 =

2.8375e+04

C2 =

1.5514e+04

C3 =

2.8375e+04

n=4时，A2为病态矩阵

n=6时：

A2=[1,1/2,1/3,1/4,1/5,1/6;1/2,1/3,1/4,1/5,1/6,1/7;1/3,1/4,1/5,1/6,1/7,1/8;1

/4,1/5,1/6,1/7,1/8,1/9;1/5,1/6,1/7,1/8,1/9,1/10;1/6,1/7,1/8,1/9,1/10,1/11] C1 =

2.9070e+07

C2 =

1.4951e+07

C3 =

2.9070e+07

北航数值分析大作业三

一、题目：

关于x, y, t, u, v, w的下列方程组

0.5cos2.670.5sin1.070.5cos3.740.5sin0.79tuvwxtuvwytuvwxtuvwy

以及关于z, t, u的下列二维数表

z u

t 0 0.4 0.8 1.2 1.6 2.0

0 -0.5 -0.34 0.14 0.94 2.06 3.5

0.2 -0.42 -0.5 -0.26 0.3 1.18 2.38

0.4 -0.18 -0.5 -0.5 -0.18 0.46 1.42

0.6 0.22 -0.34 -0.58 -0.5 -0.1 0.62

0.8 0.78 -0.02 -0.5 -0.66 -0.5 -0.02

1.0 1.5 0.46 -0.26 -0.66 -0.74 -0.5

确定了一个二元函数z = f(x, y)。

1、试用数值方法求出f(x, y)在区域 {(,)|00.8,0.51.5}Dxyxy上的一个近似表达式

,0(,)krsrsrspxycxy

要求(,)pxy一最小的k值达到以下的精度

10202700((,)(,))10ijijijfxypxy

其中，0.08,0.50.05ijxiyj。

2、计算****(,),(,)ijijfxypxy（i = 1, 2, …，8；j = 1, 2,…，5）的值，以观察(,)pxy逼近(,)fxy的效果，其中，*ix=0.1i, *jy=0.5+0.2j。

说明：1、用迭代方法求解非线性方程组时，要求近似解向量()kx满足

()(1)()12||||/||||10kkkxxx

2、作二元插值时，要使用分片二次代数插值。

3、要由程序自动确定最小的k值。

4、打印以下内容：

●算法的设计方案。

●全部源程序（要求注明主程序和每个子程序的功能）。 ●数表：,,ijxy(,)ijfxy （i = 0,1,2,…,10；j = 0,1,2,…,20）。

北航数值分析大作业一

1 / 13 《数值分析B》大作业一

SY1103120 朱舜杰

一．算法设计方案：

1.矩阵A的存储与检索

将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] .

由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素aij的方法是：

A的带内元素aij=C中的元素ci-j+2,j

2.求解λ1，λ501，λs

① 首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs；

如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。

②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求出对应的按摸最大的特征值λ，max，

如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。

3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik（k=1,2，…，39）。

使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。

4.求解A的（谱范数）条件数cond（A）2和行列式detA。

①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和最小特征值。 2 / 13 ②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。

二．源程序

#include

#define E 1.0e-12 /*定义全局变量相对误差限*/

int max2(int a,int b) /*求两个整型数最大值的子程序*/

{

if(a>b)

return a;

else

return b;

数值分析期末大作业

- 1 -

一、问题提出

设方程f(x)=x3-3x-1=0有三个实根 x*1=1.8793 ,

x*2=-0.34727 ,x*3=-1.53209现采用下面六种不同计算格式，求 f(x)=0的根 x*1 或x*2 。

1、 x = 213xx

2、 x = 313x

3、 x = 313x

4、 x = 312x

5、 x = x13

6、 x = x - 1133123xxx

二、目的和意义

1、通过实验进一步了解方程求根的算法；

2、认识选择计算格式的重要性；

3、掌握迭代算法和精度控制；

4、明确迭代收敛性与初值选取的关系。

三、结构程序设计

本程序实在matlab软件上进行操作的。首先建立一个空白的M-文件。在编辑器中输入以下内容，并保存。

function [X1,m,n,q]=shizi1(p)

x=zeros(100,1);

x=double(x); - 2 -

x(1,1)=p;

i=1;

deltax=100;

while (i<100 & deltax > 0.000001)

x(i+1,1)=(3*x(i,1)+1)/x(i,1)^2

deltax=abs(x(i+1,1)-x(i,1));

i=i+1;

end

X1=x(1,1);

m=i;

n=x(i,1);

q=deltax;

以上是运行函数，下一步在建立一个执行M-文件，输入以下内容，并保存。其中X1为初始值，m为迭代次数，n为最后得到的值，q为|xk+1-xk|。

clear all;

clc;

p=1.8;

[X1,m,n,q]=shizi1(p)

1、对第一个迭代公式，在执行文件中输入p=1.8；[X1,m,n,q]=shizi1(p)。得到如下结果如下：初值为- 3 - 1.8，迭代100次，精度为10-6。

可见该迭代公式是发散的，将初值改为-1.5，其他均条件不变。

作业分析数据

为规范动火作业安全检查分析，确保特殊作业在安全条件下有序开展，依据《化学品生产单位特殊作业安全规程》GB30871-2014第5.4.2条动火分析合格标准：当被测气体爆炸极限大于或等于4%时，其被测浓度不大于0.5%（体积分数）。进行动火前的分析要求依据下列表格中的数据进行。

一、动火作业分析数据依据：

可燃气体名称分子式爆炸极限毒性下限（LEL%）上限（LEL%）

甲烷 CH4 5 15

苯 C6H6 1.3 7.1 中毒

硫化氢 H2S 4.3 45.5 神经

一氧化碳 CO 12.5 74.2 剧毒

氢 H2 4 75

氨 NH3 16 25 低毒

注：以上数据来源于《SH3063—1999石油化工企业可燃气体和有毒气体检测报警器设计规范》，并与《常用化学危险品安全手册》进行了对照。

二、受限空间作业分析数据依据：

气体名称分子式报警值允许浓度

(ppm) 最高允许浓度

(mg/m3) 底限(ppm) 高限(ppm)

氧气 O2 ----- ---- 18%-22%

氨 NH3 39.5 47.4 30

一氧化碳 CO 24 136 50 30

苯 C6H6 2.87 6 40

硫化氢 H2S 10 15 10 10

下载提示

1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。

数值分析大作业

合集下载

北航数值分析大作业三

北航数值分析大作业一

数值分析期末大作业

作业分析数据

文档推荐

最新文档