六点格式与加权六点格式.docx
- 文档编号:27979476
- 上传时间:2023-07-07
- 格式:DOCX
- 页数:19
- 大小:315.19KB
六点格式与加权六点格式.docx
《六点格式与加权六点格式.docx》由会员分享,可在线阅读,更多相关《六点格式与加权六点格式.docx(19页珍藏版)》请在冰豆网上搜索。
六点格式与加权六点格式
偏微分方程实验三
内容:
六点格式与加权六点格式
完成时间:
2014年03月09日
格式三:
六点格式
#include
#include
#defineNk49
#defineNt10
#definer25
#defineh0.2
#definet1
voidfunc(floatA[Nk][Nk],floatU[51][11])
{
inti,k,j,l,g;
floatX[Nk][1]={0};
//将三对角矩阵A化为上三角形式
for(k=1;k<=Nt;k++)
{
g=0;
//初始化B矩阵
floatB[Nk][1]={0};
if(!
g)
{
for(l=0;l { if(l==0) { B[l][0]=r*U[l][1]/2+(1-r)*U[l+1][0]+r*U[l+2][0]/2+r*U[l][0]/2; } elseif(l==Nk-1) { B[l][0]=-r*U[l+1][Nk]/2+(1-r)*U[l][Nk-1]+r*U[l+1][Nk-1]/2+r*U[l-1][Nk-1]/2; } else B[l][0]=(1-r)*U[k][l]+r*U[k+1][l]/2+r*U[k-1][l]/2; } g=1; } //***追过程*** A[0][1]=A[0][1]/A[0][0]; B[0][0]=B[0][0]/A[0][0]; A[0][0]=1; for(i=1;i { A[i][i+1]=A[i][i+1]/(A[i][i]-U[i-1][i]*A[i][i-1]); B[i][0]=(B[i][0]-B[i-1][0]*A[i][i-1])/(A[i][i]-U[i-1][i]*A[i][i-1]); A[i][i]=1; A[i][i-1]=0; } A[Nk-1][Nk-1]=1; //***赶过程*** X[Nk-1][0]=B[Nk-1][0]; for(i=Nk-2;i>=0;i--) { X[i][0]=B[i][0]-U[i][i+1]*X[i+1][1]; } for(j=1;j<=Nk;j++) { U[j][k]=X[j-1][1]; } } for(i=0;i<51;i++) { for(j=0;j<11;j++) { printf("%13.1f",U[i][j]); } printf("\n"); } } voidmain() { inti,j,k=0; floatA[Nk][Nk]={0},U[51][11]={0},x0=0,t0=0; //三对角矩阵A的初始化 A[0][0]=1+r; A[0][1]=-r/2; A[Nk-1][Nk-2]=-r/2; A[Nk-1][Nk-1]=1+r; for(i=1;i { for(j=i-1;j<=i+1;j++) { if(i==j) A[i][j]=1+r; else A[i][j]=-r/2; } } //初始时刻的U矩阵 for(i=0;i<51;i++) U[i][0]=exp(x0+h*i); for(j=1;j<11;j++) U[0][j]=exp((t0+t)*j); for(j=1;j<11;j++) U[50][j]=exp((t0+t)*j+1); //判断A矩阵是否对角占优 for(i=0;i { floatsum=0; for(j=0;j { if(i! =j) sum+=A[i][j]; } if(A[i][i]>sum) k++; } if(k==Nk) { printf("A是严格对角占优! \n"); func(A,U); } else printf("A不是严格对角占优"); } 部分运行结果: Matlab作图: Plot(U); 加权六点格式: #include #include #defineNk49 #defineNt10 #definer25 #defineh0.2 #definet1 voidfunc(floatA[Nk][Nk],floatU[51][11],floattheta) { inti,k,j,l,g; floatX[Nk][1]={0}; //将三对角矩阵A化为上三角形式 for(k=1;k<=Nt;k++) { g=0; //初始化B矩阵 floatB[Nk][1]={0}; if(! g) { for(l=0;l { if(l==0) { B[l][0]=theta*r*U[l][1]+(1-2*(1-theta)*r)*U[l+1][0]+(1-theta)*r*U[l+2][0]+(1-theta)*r*U[l][0]; } elseif(l==Nk-1) { B[l][0]=theta*r*U[l+1][Nk]+(1-2*(1-theta)*r)*U[l][Nk-1]+(1-theta)*r*U[l+1][Nk-1]+(1-theta)*r*U[l-1][Nk-1]; } else B[l][0]=(1-2*(1-theta)*r)*U[k][l]+(1-theta)*r*U[k+1][l]+(1-theta)*r*U[k-1][l]; } g=1; } //***追过程*** A[0][1]=A[0][1]/A[0][0]; B[0][0]=B[0][0]/A[0][0]; A[0][0]=1; for(i=1;i { A[i][i+1]=A[i][i+1]/(A[i][i]-U[i-1][i]*A[i][i-1]); B[i][1]=(B[i][0]-B[i-1][0]*A[i][i-1])/(A[i][i]-U[i-1][i]*A[i][i-1]); A[i][i]=1; A[i][i-1]=0; } A[Nk-1][Nk-1]=1; //***赶过程*** X[Nk-1][0]=B[Nk-1][0]; for(i=Nk-2;i>=0;i--) { X[i][0]=B[i][0]-U[i][i+1]*X[i+1][0]; } for(j=1;j<=Nk;j++) { U[j][k]=X[j-1][1]; } } for(i=0;i<51;i++) { for(j=0;j<11;j++) { printf("%12.1f",U[i][j]); } printf("\n"); } } voidmain() { inti,j,k=0; floatA[Nk][Nk]={0},U[51][11]={0},x0=0,t0=0,theta; printf("请输入theta(0到1之间的一个数)的值: "); scanf("%f",&theta); //三对角矩阵A的初始化 A[0][0]=1+2*theta*r; A[0][1]=-r*theta; A[Nk-1][Nk-2]=-r*theta; A[Nk-1][Nk-1]=1+2*theta*r; for(i=1;i { for(j=i-1;j<=i+1;j++) { if(i==j) A[i][j]=1+2*theta*r; else A[i][j]=-r*theta; } } //初始时刻的U矩阵 for(i=0;i<51;i++) U[i][0]=exp(x0+h*i); for(j=1;j<11;j++) U[0][j]=exp((t0+t)*j); for(j=1;j<11;j++) U[50][j]=exp((t0+t)*j+1); //判断A矩阵是否对角占优 for(i=0;i { floatsum=0; for(j=0;j { if(i! =j) sum+=A[i][j]; } if(A[i][i]>sum) k++; } if(k==Nk) { printf("A是严格对角占优! \n"); func(A,U,theta); } else printf("A不是严格对角占优"); } 部分运行结果: 取theta=0.1时: Matlab作图: 取theta=0.3 作图: 取theta=0.5时为六点格式 取theta=0.8时: 作图: 取theta=0.95时: 作图:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 格式 加权