偏微分方程数值解程序作业.docx

偏微分方程数值解程序作业.docx

《偏微分方程数值解程序作业.docx》由会员分享，可在线阅读，更多相关《偏微分方程数值解程序作业.docx（9页珍藏版）》请在冰点文库上搜索。

偏微分方程数值解程序作业

偏微分方程数值解程序作业

二维热传导

#include"iostream"

#include"math.h"

#include"stdio.h"

#defineTT0.1//要求的t的值

#definePI3.1415926535

usingnamespacestd;

#defineN11

//跳点全局变量

doubledot_sheet_b[N][N];

doubledot_sheet_t[N][N];

 //显示全局变量

doubledot_sheet_old[N][N];

doubledot_sheet_new[N][N];

 doublelamda=0.5;//lamda=a*t/powf（h,2）

doubleh1=0.1,h2=0.1;

voidget_next_dot（inti,intj）;

voidinitial（void）;

voidinitial（void）;

voidinitial1（void）;

voidget_next_2_layer（void）;//每次得到两层

voidget_next_1_layer（void）;//显示一次获得一层

voidmain（void）

{

 inti,j;//行列计数使用

 initial（）;

 get_next_2_layer（）;

cout<<"output_line"<

 for（i=1;i

 {

  for（j=1;j

   printf（"%.3f|",dot_sheet_b[i][j]）;   

  cout<

   }

 printf（"**************隐式结果***************\n"）;

   initial1（）;

  get_next_1_layer（）;

  for（i=1;i

  {

   for（j=1;j

    printf（"%.3f|",dot_sheet_new[i][j]）;

   cout<

    }

 printf（"***************显示结果*******************\n"）;

}

 voidget_next_1_layer（void）//显示一次获得一层

{

 doubledelt_x2,delt_y2;

 inti,j;

 for（i=1;i

 {

  for（j=1;j

  {

   delt_x2=（dot_sheet_old[i+1][j]-2*dot_sheet_old[i][j]+dot_sheet_old[i-1][j]）;

   delt_y2=（dot_sheet_old[i][j+1]-2*dot_sheet_old[i][j]+dot_sheet_old[i][j-1]）;

      dot_sheet_new[i][j]=dot_sheet_old[i][j]+lamda*（delt_x2+delt_y2）;

  }

 }

}

voidget_next_2_layer（void）//跳点格式一次获得两层

{

 doubledelt_x2,delt_y2;

 inti,j;

  for（i=1;i

 {

  for（j=（i+1）%2+1;j

  {

   delt_x2=（dot_sheet_t[i+1][j]-2*dot_sheet_t[i][j]+dot_sheet_t[i-1][j]）;

   delt_y2=（dot_sheet_t[i][j+1]-2*dot_sheet_t[i][j]+dot_sheet_t[i][j-1]）;

   dot_sheet_b[i][j]=dot_sheet_t[i][j]+lamda*（delt_x2+delt_y2）;

//   printf（"%.3f|",dot_sheet_b[i][j]）;

  }//  cout<

 }

 for（i=1;i

 {

  for（j=（i）%2+1;j

  {

    delt_x2=（dot_sheet_t[i+1][j]-2*dot_sheet_t[i][j]+dot_sheet_t[i-1][j]）;

    delt_y2=（dot_sheet_t[i][j+1]-2*dot_sheet_t[i][j]+dot_sheet_t[i][j-1]）;

    dot_sheet_b[i][j]=lamda*（delt_x2+delt_y2）+dot_sheet_t[i][j];

//    printf（"%.3f|",dot_sheet_b[i][j]）;

  }

//  cout<

 }

  cout<<"***********************************************************"<

 for（i=1;i

 {

  for（j=1;j

  {

   dot_sheet_t[i][j]=2*dot_sheet_b[i][j]-dot_sheet_t[i][j];

//   printf（"%.4f|",dot_sheet_t[i][j]）;

   //第二层输出屏蔽，暂时只输出第一层

  }

//  cout<

 }

}

voidinitial（void）//跳点格式初始化

{

 for（inti=0;i

  dot_sheet_t[i][0]=dot_sheet_t[i][N-1]=0;

 for（intj=0;j

  dot_sheet_t[0][j]=dot_sheet_t[N-1][j]=0;

 for（i=1;i

 {

  for（j=1;j

  {

   dot_sheet_b[i][j]=dot_sheet_t[i][j]=sin（（h1*i+h2*j）*PI）;

   printf（"%.4f|",dot_sheet_t[i][j]）;

  }

   cout<

 }

 printf（"************初始化底层**********************\n"）;

}

 voidinitial1（void）//显示初始化

{

 for（inti=0;i

  dot_sheet_old[i][0]=dot_sheet_old[i][N-1]=0;

 for（intj=0;j

  dot_sheet_old[0][j]=dot_sheet_old[N-1][j]=0;

 for（i=1;i

 {

  for（j=1;j

  {

   dot_sheet_old[i][j]=sin（（h1*i+h2*j）*PI）;

  }

 }

 printf（"**************初始化底层1*******************\n"）;

}

一维热传导-----向前差分

#include"iostream"

#include"math.h"

#include"stdio.h"

#defineTT0.1//要求的t的值

#definePI3.1415926535

usingnamespacestd;

voidmain（void）

{

 doubler=0.1;//0.5

 doublestep_h=0.1;

 doublestep_t;//初始化无任何意义

 introw1,cln1,row_max;

 step_t=r*（step_h*step_h）;

// printf（"row_max=%f\n",step_t）;

来源：

（-向前差分格式_爱学_新浪博客 row_max=TT/step_t+1;//认为是第rowmax行

 printf（"row_max=%d\n",row_max）;

 doublematrix_u[20][11];//一行共11个点

//以下为带入边界条件

 for（cln1=0;cln1<=10;cln1++）

  matrix_u[0][cln1]=sin（PI*step_h*cln1）;

 for（row1=0;row1<=row_max;row1++）

  matrix_u[row1][0]=matrix_u[row1][10]=0;

 for（row1=1;row1<=row_max;row1++）

  for（cln1=1;cln1<=9;cln1++）

   matrix_u[row1][cln1]=r*（matrix_u[row1+1][cln1-1]+matrix_u[row1-1][cln1-1]）+（1-r）*（matrix_u[row1][cln1-1]）;

   for（cln1=0;cln1<=10;cln1++）

 {

  printf（"matrix_u[%d][%d]=%.5f\n",row1-1,cln1,matrix_u[row1-1][cln1]）;

  printf（"theexactvalue=%.5f\n",exp（-PI*PI*（row_max-1）*step_t）*sin（PI*step_h*cln1））;

  }

}

向后差分格式

#include"iostream"

#include"math.h"

#include"stdio.h"

#include"malloc.h"

#defineTT0.1//要求的t的值

#definePI3.1415926535

usingnamespacestd;

voidzhuigan（double*A,double*B,double*C,double*D,intn）//n表示矩阵的维数，就是行数。

{

 inti;

 B[1-1]=B[1-1];

 D[1-1]=D[1-1]/B[1-1];

  for（i=2;i<=n;i++）

  {

   C[i-2]=C[i-2]/B[i-2];

   B[i-1]=B[i-1]-A[i-1]*C[i-2];

   D[i-1]=（D[i-1]-A[i-1]*D[i-2]）/B[i-1];

来源：

（-偏微分向后差分格式_爱学_新浪博客 }

  for（i=n-1;i>=1;i--）

   D[i-1]-=C[i-1]*D[i];

}//endoffunction

voidmain（void）

{

 doubler=0.1;//0.5

 doublestep_h=0.1;

 doublestep_t;//

 introw1,cln1,row_max;

 step_t=r*（step_h*step_h）;

 printf（"step_t=%f\n",step_t）;

 row_max=TT/step_t+1;//认为是第rowmax行

 printf（"row_max=%d\n",row_max）;

//追赶系数矩阵赋值

 doublem_a[9],m_b[9],m_c[9],m_d[9];

 doublem_a_temp[9],m_b_temp[9],m_c_temp[9];//必须的

// doublel_l[9],r_r[9],r_y[9];

 for（row1=1;row1<=9;row1++）

 {m_b[row1-1]=1+2*r;}

 for（row1=2;row1<=9;row1++）

  m_a[row1-1]=-r;

 for（row1=1;row1<=9-1;row1++）

  m_c[row1-1]=-r;

 for（cln1=1;cln1<=9;cln1++）

  m_d[cln1-1]=sin（PI*step_h*cln1）;

 for（row1=1;row1

 { 

  for（cln1=1;cln1<=9;cln1++）

  {

   m_a_temp[cln1-1]=m_a[cln1-1];

   m_b_temp[cln1-1]=m_b[cln1-1];

   m_c_temp[cln1-1]=m_c[cln1-1];

  }

  zhuigan（m_a_temp,m_b_temp,m_c_temp,m_d,9）;

 }

 for（row1=1;row1<=9;row1++）

 {

  printf（"d%.4f",m_d[row1-1]）;

 }

 cout<

  cout<<"theexactvalue="<

 for（cln1=1;cln1<=9;cln1++）

  printf（"E%.4f",exp（-1*PI*PI*（row_max-1）*step_t）*sin（PI*step_h*cln1））;

}//endofmain

六点格式

#include"iostream"

#include"math.h"

#include"stdio.h"

#include"malloc.h"

#defineTT0.1//要求的t的值

#definePI3.1415926535

usingnamespacestd;

voidzhuigan（double*A,double*B,double*C,double*D,intn）

{

 inti;

 B[1-1]=B[1-1];

  D[1-1]=D[1-1]/B[1-1];

 for（i=2;i<=n;i++）

 {

  C[i-2]=C[i-2]/B[i-2];

    B[i-1]=B[i-1]-A[i-1]*C[i-2];

    D[i-1]=（D[i-1]-A[i-1]*D[i-2]）/B[i-1];

   }

 for（i=n-1;i>=1;i--）

  D[i-1]-=C[i-1]*D[i];

}//endoffunction

来源：

（-偏微分方程六点法 crank-nicolson_爱学_新浪博客voidmain（void）

{

 doubler=0.1;//0.5

 doublestep_h=0.1;

 doublestep_t;//

 introw1,cln1,row_max;

 step_t=r*（step_h*step_h）;

 printf（"step_t=%f\n",step_t）;

 row_max=TT/step_t+1;//认为是第rowmax行

 printf（"row_max=%d\n",row_max）;

  //追赶系数矩阵赋值

 doublem_a[9],m_b[9],m_c[9],m_d[9],u_d[11];//u_d是要求的值，m_d使方程用的

 doublem_a_temp[9],m_b_temp[9],m_c_temp[9];//必须的

  // doublel_l[9],r_r[9],r_y[9];

 for（row1=1;row1<=9;row1++）

  m_b[row1-1]=1+r;

 for（row1=2;row1<=9;row1++）

  m_a[row1-1]=-r/2;

 for（row1=1;row1<=9-1;row1++）

  m_c[row1-1]=-r/2;

 u_d[0]=u_d[10]=0;

 for（cln1=1;cln1<=9;cln1++）

  u_d[cln1]=sin（PI*step_h*cln1）;//初始化

  for（row1=1;row1<=row_max;row1++）//第100行为结果

 {

  for（cln1=1;cln1<=9;cln1++）

   m_d[cln1-1]=r/2*（u_d[cln1-1]）+（1-r）*u_d[cln1]+r/2*（u_d[cln1+1]）;

  for（cln1=1;cln1<=9;cln1++）

  {

   m_a_temp[cln1-1]=m_a[cln1-1];

   m_b_temp[cln1-1]=m_b[cln1-1];

   m_c_temp[cln1-1]=m_c[cln1-1];

  }

  zhuigan（m_a_temp,m_b_temp,m_c_temp,m_d,9）;

  for（cln1=1;cln1<=9;cln1++）

   u_d[cln1]=m_d[cln1-1];

 }

  for（row1=1;row1<=9;row1++）

 {

  printf（"d%.4f",m_d[row1-1]）;

 }

  cout<

  cout<<"theexactvalue="<

  for（cln1=1;cln1<=9;cln1++）

  printf（"E%.4f",exp（-1*PI*PI*（99）*step_t）*sin（PI*step_h*cln1））;

 }//endofmain

来源：

（-偏微分方程六点法 crank-nicolson_爱学_新浪博客来源：

（-偏微分向后差分格式_爱学_新浪博客

来源：

（-向前差分格式_爱学_新浪博客来源：

（-二维热传导的c程序实现方法，偏微分作业_爱学_新浪博客