数值分析实验报告.docx
《数值分析实验报告.docx》由会员分享，可在线阅读，更多相关《数值分析实验报告.docx（18页珍藏版）》请在冰点文库上搜索。
数值分析实验报告
数值分析报告
班级：
计算机14XX班
学号:
201437XX
姓名：
XX
（1）
其中
为了使种群中各个年龄的种群在数量上维持不变则：
（2）由题可得矩阵
解方程得到
x1=8481.012658x2=2892.405063x3=1335.443038x4=601.265823x5=140.506329
（3）将H修改为H=（500,500,500,500,500）（注：
将代码中的H[5]修改为H[5]={500,500,500,5000,500}即可）
得到：
x1=+11772.151899x2=+4208.860759x3=+2025.316456x4=+715.189873x5=-213.924051
x5为负值，所以在
（2）问的条件下无法满足条件
此时需要修改出生率或生存率来达到
1．改变各年龄阶段物种的生存率
此时需要减少高年龄段的死亡率，才能够使年龄大的物种数量变大。
如果我们只改变s4的值，一直改到1，计算的结果为：
x1=11772.151899x2=4208.860759x3=2025.316456x4=+715.189873x5=215.189873
由此可见x5<500不符合题意
于是同时增加S4与S5的值
直到S[0.3,0.6,0.8,0.8]
得到：
x1=23855.421687x2=6656.626506x3=3493.975904x4=2295.180723x5=1336.144578
从中可以看出X5>500，符合题意。
但是从生物学的角度来看，低年龄段的种群数量过于多，容易受到天敌或是自然灾害的袭击不符合现实。
2．改变出生率
改变出生率b3=4,b4=2得到
x1=22822.580645x2=8629.032258x3=4677.419355x4=2306.451613x5=422.580645
X5<500不符合题意
使b3=4,b4=1得到：
x1=45000.000000x2=17500.000000x3=+10000.000000x4=5500.000000x5=1700.000000
X5>500符合题意
比起第一种情况各个年龄段的数量是相对平均的，但是各个年龄段的种群数量很大，尤其是第一年龄段，从中可以看出一个种群想要蓬勃的发展小年龄段的种群数量就必须要多。
附页：
1.LU分解求逆矩阵算法
矩阵的一种有效而广泛应用的分解方法是矩阵的LU三角分解，将一个n阶矩阵A分解为一个下三角矩阵L和一个上三角矩阵U的乘积。
所以首先对矩阵进行三角分解，这里采用Doolittle分解，即分解为一个下三角矩阵（对角元素为1），和一个上三角矩阵的乘积。
再进行相应的处理。
所以，矩阵求逆的算法流程可表述如下：
图1矩阵求逆流程图
1）进行LU分解；
2）对分解后的L阵（下三角矩阵）和U阵（上三角矩阵）进行求逆;；
3）L阵的逆矩阵和U阵的逆矩阵相乘，即可求得原来矩阵的逆。
即：
（1）
1.1矩阵的LU分解
若n阶方阵
的各阶顺序主子式不等于零，即：
（2）
则A的LU分解
存在且唯一。
（3）
由矩阵的乘法原理,可推导出LU分解的迭代算法
（4）
（5）
（6）
（7）
矩阵的LU分解是一个循环迭代的过程,U矩阵是从第1行迭代到第n行,而L矩阵则是从第1列迭代到第n列,且U矩阵先于L矩阵一个节拍。
1.2L矩阵和U矩阵求逆
首先假设下三角矩阵L的逆矩阵为
，不失一般性，考虑4阶的情况，利用
，有：
（1）
，
，
；
（2）
（3）
（4）
。
从而求得下三角矩阵L的逆矩阵R式如下：
，（8）
上三角矩阵U的逆矩阵可以由下式得到：
。
，（9）
矩阵求逆是一个迭代的过程，依次循环,迭代
次,求出整个逆矩阵。
其中U矩阵的循环迭代时按行顺序，列倒序进行，L矩阵的循环迭代按列顺序，行顺序进行，直到计算出整个矩阵的所有结果为止。
2.求解线性方程的源代码
#include
#include
#include
//采用LU分解求矩阵的逆
doublea[5][5]={0};//定义矩阵A-E的逆
doubleL[5][5]={1,0,0,0,0,
0,1,0,0,0,
0,0,1,0,0,
0,0,0,1,0,
0,0,0,0,1};//定义L矩阵
doubleU[5][5]={0};//定义U矩阵
doubleA[5][5]={-1,0.4,0,0,0,
0,-1,0.6,0,0,
5,0,-1,0.6,0,
3,0,0,-1,0.4,
0,0,0,0,-1};//定义矩阵A-E
doubleH[5]={0,500,400,200,100};//定义矩阵H=[0，b1,b2,b3,b4]
double_L[5][5]={0};//定义L矩阵的逆矩阵
double_U[5][5]={0};//定义U矩阵的逆矩阵
doubleX[5]={0};//定义矩阵X
//
voidcalLU（）//计算LU矩阵
{
inti,j,k,t;
for（i=0;i<5;i++）
for（j=0;j<5;j++）
{
if（i==0）
U[i][j]=A[i][j];
else
{
for（t=0;t{
if（t==0）
L[i][t]=A[i][t]/U[t][t];
if（t==1）
L[i][t]=（A[i][t]-L[i][0]*U[0][t]）/U[t][t];
if（t==2）
L[i][t]=（A[i][t]-L[i][0]*U[0][t]-L[i][1]*U[1][t]）/U[t][t];
if（t==3）
L[i][t]=（A[i][t]-L[i][0]*U[0][t]-L[i][1]*U[1][t]-L[i][2]*U[2][t]）/U[t][t];
if（t==4）
L[i][t]=（A[i][t]-L[i][0]*U[0][t]-L[i][1]*U[1][t]-L[i][2]*U[2][t]-L[i][3]*U[3][t]）/U[t][t];
}
//计算U矩阵
if（i==1）
{
for（k=1;k<=5-i;k++）
U[1][k]=A[1][k]-L[1][0]*U[0][k];
}
if（i==2）
{
for（k=2;k<=6-i;k++）
U[2][k]=A[2][k]-L[2][0]*U[0][k]-L[2][1]*U[1][k];
}
if（i==3）
{
for（k=3;k<=7-i;k++）
U[3][k]=A[3][k]-L[3][0]*U[0][k]-L[3][1]*U[1][k]-L[3][2]*U[2][k];
}
if（i==4）
{
for（k=4;k<=8-i;k++）
U[4][k]=A[4][k]-L[4][0]*U[0][k]-L[4][1]*U[1][k]-L[4][2]*U[2][k]-L[4][3]*U[3][k];
}
}
}
}
voidoutputLU（）//输出LU矩阵
{
printf（"L矩阵:
"）;
for（inti=0;i<5;i++）
{
printf（"\n"）;
for（intj=0;j<5;j++）
{
printf（"%+.3lf",L[i][j]）;
}
}
printf（"\nU矩阵:
"）;
for（intk=0;k<5;k++）
{
printf（"\n"）;
for（intn=0;n<5;n++）
{
printf（"%+.3lf",U[k][n]）;
}
}
printf（"\n"）;
}
void_calLU（）//计算LU的逆矩阵
{
inti,j,k;
for（i=0;i<5;i++）
for（j=0;j<5;j++）
{
if（i==j）
{
_L[i][j]=1/L[i][j];
_U[i][j]=1/U[i][j];
}
}
for（i=0;i<5;i++）//求L的逆
for（j=0;j<5;j++）
{
if（i>j）
{
_L[j][i]=0;
}
elseif（i{
_U[j][i]=0;
doublet=L[j][j],m=0;
for（k=i;km=m+_L[k][i]*L[j][k];
_L[j][i]=-1/t*m;
}
}
for（i=4;i>=0;i--）//计算U逆矩阵
for（j=4;j>=0;j--）
{
if（i>j）
{
doublet=U[j][j],m=0;
for（k=j+1;k<=i;k++）
m=m+U[j][k]*_U[k][i];
_U[j][i]=-1/t*m;
}
}
}
void_outputLU（）//输出LU逆矩阵
{
printf（"L逆矩阵:
"）;
for（inti=0;i<5;i++）
{
printf（"\n"）;
for（intj=0;j<5;j++）
{
printf（"%+.3lf",_L[i][j]）;
}
}
printf（"\nU逆矩阵:
"）;
for（intk=0;k<5;k++）
{
printf（"\n"）;
for（intn=0;n<5;n++）
{
printf（"%+.3lf",_U[k][n]）;
}
}
}
voidcalA（）//计算A的逆矩阵A=_U*_L
{
inti,j,m;
for（i=0;i<5;i++）
for（j=0;j<5;j++）
{
for（m=0;m<5;m++）
a[i][j]=a[i][j]+_U[i][m]*_L[m][j];
}
}
voidoutputcalA（）//输出逆矩阵
{printf（"\n逆矩阵:
"）;
for（inti=0;i<5;i++）
{
printf（"\n"）;
for（intj=0;j<5;j++）
{
printf（"%+.3lf",a[i][j]）;
}
}
}
voidcal（）//计算结
{
inti,j,k;
for（i=0;i<5;i++）
{
for（k=0;k<5;k++）
X[i]=X[i]+H[k]*a[k][i];
}
}
voidoutputcal（）//输出结果
{
printf（"\n结果:
"）;
for（intj=0;j<5;j++）
{
printf（"x%d=%+.6lf",j+1,X[j]）;
}
printf（"\n"）;
}
intmain（）
{
calLU（）;//计算LU矩阵
outputLU（）;//输出LU矩阵
_calLU（）;//计算LU矩阵的逆矩阵
_outputLU（）;//输出LU逆矩阵
calA（）;//计算逆矩阵
outputcalA（）;//输出逆矩阵
cal（）;//计算结果
outputcal（）;//输出结果
return0;
}
升级会员 