北京理工大学徐特立学院数值分析大作业上机实验.docx
- 文档编号:25725794
- 上传时间:2023-06-11
- 格式:DOCX
- 页数:35
- 大小:228.57KB
北京理工大学徐特立学院数值分析大作业上机实验.docx
《北京理工大学徐特立学院数值分析大作业上机实验.docx》由会员分享,可在线阅读,更多相关《北京理工大学徐特立学院数值分析大作业上机实验.docx(35页珍藏版)》请在冰豆网上搜索。
北京理工大学徐特立学院数值分析大作业上机实验
北京理工大学徐特立学院数值分析课后上机实验选做
教材:
数值计算方法(2011第一版).丁丽娟,程杞元.高等教育出版社
^以下代码作者原创^
超链接:
1.2(题目,原理,截图,代码)
2.2(题目,原理,截图,代码)
3.1(题目,原理,截图,c代码,m代码)
5.1(题目,原理,截图,代码)
5.3(题目,原理,截图,代码)
第一章:
数值计算中的误差
2、题目简介:
利用pi/4=1-1/3+1/5-1/7。
。
。
级数计算pi的近似值。
输入:
误差值
输出:
求和项数,并输出pi值
工具:
C语言
运行环境:
VC-6.0
计算公式及原理:
利用pi/4=1-1/3+1/5-1/7。
。
。
级数计算pi的近似值,由数学原理可知误差会小于首次舍弃的项,可以编写循环实现。
程序运行结果截图:
程序代码:
(c语言)
#include
voidmain()
{
printf("第一章第2题求pi,欢迎使用,请按提示操作。
\n");
inti=1,n=0,k=1;
doublee,pi,er;
printf("请输入误差(例如1e-4):
");
scanf("%lf",&e);
printf("请稍候。
。
。
\n");
er=e;
pi=0;
while(er>=e)
{
pi+=k*1.0/i;k=-k;er=1.0/i;
i+=2;n++;
}
pi*=4;
printf("%d项求和后可以达到%.10lf精度,这时pi=%.10lf\n",n,e,pi);
getchar();
getchar();
}
第二章:
解线性方程组的直接方法
2、题目简介:
用MATLAB软件编程实现追赶法求解三对角方程组的算法,并考虑梯形电阻电路问题,电路如下:
其中电路中的各个电流{
,
,…,
}须满足下列线性方程组:
设
,
,运用求各段电路的电流量。
输入:
三列对角线元素向量,右侧常数元素向量
输出:
三对角方程的解
工具:
m语言
运行环境:
MATLABR2012.b
计算公式与原理:
上述方程组可用矩阵表示为:
根据三对角方程追赶法原理以及公式可解。
程序运行结果截图:
程序代码:
(matlab)
%说明:
追赶法是适用于三对角矩阵的线性方程组求解的方法,并不适用于其他类型矩阵。
%定义三对角矩阵A的各组成单元。
方程为Ax=d
functionx=zhuiganfa
a=input('输入左下角对角元素a(n)(例如[123]):
');
b=input('输入对角元素b(n)(例如[1234]):
');
c=input('输入左上角对角元素c(n)(例如[123]):
');
d=input('输入列矩阵d(n)例如[1234]:
');
n=length(b);
u0=0;y0=0;
%“追”的过程
L
(1)=b
(1);
y
(1)=(d
(1))/L
(1);
u
(1)=c
(1)/L
(1);
fori=2:
(n-1)
L(i)=b(i)-a(i-1)*u(i-1);
y(i)=(d(i)-y(i-1)*a(i-1))/L(i);
u(i)=c(i)/L(i);
end
L(n)=b(n)-a(n-1)*u(n-1);
y(n)=(d(n)-y(n-1)*a(n-1))/L(n);
%“赶”的过程
x(n)=y(n);
fori=(n-1):
-1:
1
x(i)=y(i)-u(i)*x(i+1);
end
第三章:
解线性方程组的迭代法
1、题目简介:
试分别用
(1)Jacobi迭代法;
(2)Gauss-Seidel迭代法;(3)共轭梯度法解线性方程组
迭代初始向量取
=
0,0,0,0,0
。
输入:
根据要求输入精度或者步数,系数矩阵,向量b,初值向量等
输出:
解的精度,方程的解,要达到精度所需步数
工具:
C语言/m语言
运行环境:
VC-6.0/MATLAB
计算公式及原理:
根据相关迭代原理和算法,设计程序求解,程序可以按步数测试,也可以按精度测试。
运行结果截图:
(C语言)
运行结果截图(MATLAB)
程序代码:
(c语言)
Jacobi迭代法解线性方程组:
#include
#include
voidfun(doublea[][100],double*b,double*x,double*y,double*max,intN)
{
inti,j;
doublesum;
for(i=1;i<=N;i++)
x[i]=y[i];
for(i=1;i<=N;i++)
{
sum=0;
for(j=1;j<=N;j++)
{
if(j!
=i)
sum+=(a[i][j])*(x[j]);
}
y[i]=(b[i]-sum)/a[i][i];
}
*max=fabs(y[1]-x[1]);
for(i=1;i<=N;i++)
if(fabs(y[i]-x[i])>*max)
*max=fabs(y[i]-x[i]);
}
voidmain()
{
printf("欢迎来到Jacobi迭代法解线性方程组,请先确保该法收敛,按提示操作*^_^*\n");
intN,i,j,k=0,mode,num;
doublea[100][100],b[100],x[100],y[100],e,sum,max;
printf("请输入结束条件,迭代次数限制请按1,误差限制请按2:
");
scanf("%d",&mode);
if(mode==2)
{
printf("请输入误差限e:
");
scanf("%lf",&e);
}
else
{
printf("请输入迭代次数:
");
scanf("%d",&num);
}
printf("请输入方程阶数:
");
scanf("%d",&N);
printf("请输入系数矩阵:
\n");
for(i=1;i<=N;i++)
for(j=1;j<=N;j++)
scanf("%lf",&a[i][j]);
printf("请输入向量b:
\n");
for(i=1;i<=N;i++)
scanf("%lf",&b[i]);
printf("请输入初始解x:
");
for(j=1;j<=N;j++)
scanf("%lf",&x[j]);
for(i=1;i<=N;i++)
{
sum=0;
for(j=1;j<=N;j++)
{
if(j!
=i)
sum+=a[i][j]*x[j];
}
y[i]=(b[i]-sum)/a[i][i];
}
max=fabs(y[1]-x[1]);
for(i=1;i<=N;i++)
if(fabs(y[i]-x[i])>max)
max=fabs(y[i]-x[i]);
printf("%d",k);
for(i=1;i<=N;i++)
printf("x[%d]=%-10.7lf",i,x[i]);
printf("\n");
k++;
printf("%d",k);
for(i=1;i<=N;i++)
printf("x[%d]=%-10.7lf",i,y[i]);
printf("\n");
if(mode==1)
{
while(k { fun(a,b,x,y,&max,N); k++; if(k<10) printf("%d",k); else printf("%d",k); for(i=1;i<=N;i++) printf("x[%d]=%-10.7lf",i,y[i]); printf("\n"); } } else { while(max>e) { fun(a,b,x,y,&max,N); k++; if(k<10) printf("%d",k); else printf("%d",k); for(i=1;i<=N;i++) printf("x[%d]=%-10.7lf",i,y[i]); printf("\n"); } } getchar(); getchar(); } Gauss-Sdidel迭代法解线性方程组: #include #include voidfun(doublea[][100],double*b,double*x,double*y,double*max,intN) { inti,j; doublesum,ee[100]; for(i=1;i<=N;i++) x[i]=y[i]; for(i=1;i<=N;i++) { sum=0; for(j=1;j<=N;j++) { if(j! =i) sum+=(a[i][j])*(x[j]); } y[i]=(b[i]-sum)/a[i][i]; ee[i]=y[i]-x[i]; x[i]=y[i]; } *max=fabs(ee[1]); for(i=1;i<=N;i++) if(fabs(ee[i])>*max) *max=fabs(ee[i]); } voidmain() { printf("欢迎来到Gauss-Sdidel迭代法解线性方程组,请先确保该法收敛,按提示操作*^_^*\n"); intN,i,j,k=0,mode,num; doublea[100][100],b[100],x0[100],x[100],y[100],e,ee[100],sum,max; printf("请输入结束条件,迭代次数限制请按1,误差限制请按2: "); scanf("%d",&mode); if(mode==2) { printf("请输入误差限e: "); scanf("%lf",&e); } else { printf("请输入迭代次数: "); scanf("%d",&num); } printf("请输入方程阶数: "); scanf("%d",&N); printf("请输入系数矩阵: \n"); for(i=1;i<=N;i++) for(j=1;j<=N;j++) scanf("%lf",&a[i][j]); printf("请输入向量b: \n"); for(i=1;i<=N;i++) scanf("%lf",&b[i]); printf("请输入初始解x: "); for(j=1;j<=N;j++) { scanf("%lf",&x[j]); x0[j]=x[j]; } for(i=1;i<=N;i++) { sum=0; for(j=1;j<=N;j++) { if(j! =i) sum+=a[i][j]*x[j]; } y[i]=(b[i]-sum)/a[i][i]; ee[i]=y[i]-x[i]; x[i]=y[i]; } max=fabs(ee[1]); for(i=1;i<=N;i++) if(fabs(ee[i])>max) max=fabs(ee[i]); printf("%d",k); for(i=1;i<=N;i++) printf("x[%d]=%-10.7lf",i,x0[i]); printf("\n"); k++; printf("%d",k); for(i=1;i<=N;i++) printf("x[%d]=%-10.7lf",i,x[i]); printf("\n"); if(mode==1) { while(k { fun(a,b,x,y,&max,N); k++; if(k<10) printf("%d",k); else printf("%d",k); for(i=1;i<=N;i++) printf("x[%d]=%-10.7lf",i,y[i]); printf("\n"); } } else { while(max>e) { fun(a,b,x,y,&max,N); k++; if(k<10) printf("%d",k); else printf("%d",k); for(i=1;i<=N;i++) printf("x[%d]=%-10.7lf",i,y[i]); printf("\n"); } } getchar(); getchar(); } 共轭梯度法c语言源程序: #include #include #defineN100 voidfun1(doublea[][N],double*b,double*c) { inti,j; doublesum=0; for(i=1;i { sum=0; for(j=1;j sum+=a[i][j]*b[j]; c[i]=sum; } } doublefun2(double*c,double*d) { inti; doublesum=0; for(i=1;i sum+=c[i]*d[i]; returnsum; } voidmain() { inti,j,mode,n,num,k=0; doublee,a[N][N]={0}; doubleb[N]={0},x0[N]={0},x[N]={0},r[N]={0},d[N]={0}; doubleBeta,Lambda,max=0,c1[N]={0},c2[N]={0}; printf("欢迎来到共轭梯度法解线性方程组,请先确保该法收敛,按提示操作*^_^*\n"); printf("请输入结束条件,迭代次数限制请按1,误差限制请按2: "); scanf("%d",&mode); if(mode==2) { printf("请输入误差限e: "); scanf("%lf",&e); } else { printf("请输入迭代次数: "); scanf("%d",&num); } printf("请输入方程阶数: "); scanf("%d",&n); printf("请输入系数矩阵: \n"); for(i=1;i<=n;i++) for(j=1;j<=n;j++) scanf("%lf",&a[i][j]); printf("请输入向量b: \n"); for(i=1;i<=n;i++) scanf("%lf",&b[i]); printf("请输入初始解x: "); for(j=1;j<=n;j++) { scanf("%lf",&x[j]); x0[j]=x[j]; } fun1(a,x,c1); for(i=1;i<=n;i++) { d[i]=b[i]-c1[i]; r[i]=d[i]; max=r[1]; if(fabs(r[i])>max) max=fabs(r[i]); } fun1(a,d,c1); Lambda=fun2(r,d)/fun2(d,c1); printf("%d",k); for(i=1;i<=n;i++) printf("x[%d]=%-10.7lf",i,x0[i]); printf("\n"); k++; printf("%d",k); for(i=1;i<=n;i++) { x[i]=x[i]+Lambda*d[i]; printf("x[%d]=%-10.7lf",i,x[i]); } printf("\n"); if(mode==1) { while(k { k++; fun1(a,x,c1); for(i=1;i<=n;i++) { r[i]=b[i]-c1[i]; max=r[1]; if(fabs(r[i])>max) max=fabs(r[i]); } fun1(a,d,c2); Beta=-fun2(r,c2)/fun2(d,c2); for(i=1;i<=n;i++) d[i]=r[i]+Beta*d[i]; fun1(a,d,c2); Lambda=fun2(r,d)/fun2(d,c2); if(k<10) printf("%d",k); else printf("%d",k); for(i=1;i<=n;i++) { x[i]+=Lambda*d[i]; printf("x[%d]=%-10.7lf",i,x[i]); } printf("\n"); } } else { while(max>e) { k++; fun1(a,x,c1); for(i=1;i<=n;i++) { r[i]=b[i]-c1[i]; max=r[1]; if(fabs(r[i])>max) max=fabs(r[i]); } fun1(a,d,c2); Beta=-fun2(r,c2)/fun2(d,c2); for(i=1;i<=n;i++) d[i]=r[i]+Beta*d[i]; fun1(a,d,c2); Lambda=fun2(r,d)/fun2(d,c2); if(k<10) printf("%d",k); else printf("%d",k); for(i=1;i<=n;i++) { x[i]+=Lambda*d[i]; printf("x[%d]=%-10.7lf",i,x[i]); } printf("\n"); } } getchar(); getchar(); } 程序代码(matlab) (1)Jacobi迭代法MatLab源程序。 A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15]; b=[12,-27,14,-17,12]; x0=[0,0,0,0,0];x1=x0; k=0; fori=1: 5 sum=0; forj=1: 5 ifj~=i sum=sum+A(i,j)*x0(j); end; end; x1(i)=(b(i)-sum)/A(i,i); end; whileabs(norm(x1-x0,inf))>1e-8&k x0=x1; fori=1: 5 sum=0; forj=1: 5 ifj~=i sum=sum+A(i,j)*x0(j); end; end; x1(i)=(b(i)-sum)/A(i,i); end; k=k+1; end; fprintf('精度'); 1e-8 fprintf('解'); x1 fprintf('迭代次数'); k 输出结果: x1= 1.000000024731627-2.0000000226631672.999999958446139-1.9999999933358230.999999970325168 k=88 迭代了88次。 (2)Gauss-Seidel迭代法MatLab源程序。 A=[10,1,2,3,4;1,9,-1,2,-3;2,-1,7,3,-5;3,2,3,12,-1;4,-3,-5,-1,15];b=[12,-27,14,-17,12]; x0=[0,0,0,0,0];x1=x0; Nmax=1000; k=1;sum=0; forj=2: 5 sum=sum+A(1,j)*x0(j); end; x1 (1)=(b (1)-sum)/A(1,1); fori=2: 4 sum=0; forj=1: (i-1) sum=sum+A(i,j)*x0(j); end; x1(i)=b(i)-sum; sum=0; forj=(i+1): 5 sum=sum+A(i,j)*x0(j); end; x1(i)=(x1(i)-sum)/A(i,i); end; sum=0; forj=1: 4 sum=sum+A(5,j)*x0(j); end; x1(5)=(b(5)-sum)/A(5,5); whileabs(norm(x1-x0,inf))>1e-8&k x0=x1; sum=0; forj=2: 5 sum=sum+A(1,j)*x0(j); end; x1 (1)=(b (1)-sum)/A(1,1); fori=2: 4 sum=0; forj=1: (i-1) sum=sum+A(i,j)*x0(j); end; x1(i)=b(i)-sum; sum=0; forj=(i+1): 5 sum=sum+A(i,j)*x0(j); end; x1(i)=(x1(i)-sum)/A(i,i); end; sum=0; forj=1: 4 sum=sum+A(5,j)*x0(j); end; x1(5)=(b(5)-sum)/A(5,5); k=k+1; end; fprintf('精度'); 1e-8 fprintf('解'); x1 fprintf('迭代次数'); k 输出结果: x1= 1.000000024731627-2.0000000226631662.999999958446139
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 北京理工大学 徐特立 学院 数值 分析 作业 上机 实验