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

二维热传导

#include "iostream"

#include "math.h"

#include "stdio.h"

#define TT 0.1//要求的t的值

#define PI 3.1415926535

using namespace std;

#define N 11

//跳点全局变量

double dot_sheet_b[N][N];

double dot_sheet_t[N][N];

//显示全局变量

double dot_sheet_old[N][N];

double dot_sheet_new[N][N];

double lamda=0.5;//lamda=a*t/powf(h,2) double h1=0.1,h2=0.1;

void get_next_dot(int i,int j);

void initial(void);

void initial1(void);

void get_next_2_layer(void);//每次得到两层void get_next_1_layer(void);//显示一次获得一层

void main(void)

{

int i,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"); }

void get_next_1_layer(void)//显示一次获得一层

{

double delt_x2,delt_y2;

int i,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);

}

void get_next_2_layer(void)//跳点格式一次获得两层

{

double delt_x2,delt_y2;

int i,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];

}

// 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<

}

void initial(void)//跳点格式初始化

{

for(int i=0;i

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

for (int j=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);

}

cout<

}

printf("************初始化底层**********************\n"); }

void initial1(void)//显示初始化

{

for(int i=0;i

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

for (int j=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"

#define TT 0.1//要求的t的值

#define PI 3.1415926535

using namespace std;

void main(void)

{

double r=0.1;//0.5

double step_h=0.1;

double step_t;//初始化无任何意义

int row1,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);

double matrix_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("the exact value=%.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"

#define TT 0.1//要求的t的值

#define PI 3.1415926535

using namespace std;

void zhuigan(double *A,double *B,double *C,double *D,int n)//n表示矩阵的维数，就是行数。{

int i;

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];

}//end of function

void main(void)

{

double r=0.1;//0.5

double step_h=0.1;

double step_t;//

int row1,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);

//追赶系数矩阵赋值

double m_a[9],m_b[9],m_c[9],m_d[9];

double m_a_temp[9],m_b_temp[9],m_c_temp[9];//必须的

// double l_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<<"the exact value="<

for(cln1=1;cln1<=9;cln1++)

printf("E%.4f ",exp(-1*PI*PI*(row_max-1)*step_t)*sin(PI*step_h*cln1)); }//end of main

六点格式

#include "iostream"

#include "math.h"

#include "stdio.h"

#include "malloc.h"

#define TT 0.1//要求的t的值

#define PI 3.1415926535

using namespace std;

void zhuigan(double *A,double *B,double *C,double *D,int n)

{

int i;

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];

}//end of function

void main(void)

{

double r=0.1;//0.5

double step_h=0.1;

double step_t;//

int row1,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);

//追赶系数矩阵赋值

double m_a[9],m_b[9],m_c[9],m_d[9],u_d[11];//u_d 是要求的值，m_d使方程用的 double m_a_temp[9],m_b_temp[9],m_c_temp[9];//必须的

// double l_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<<"the exact value="<

for(cln1=1;cln1<=9;cln1++)

printf("E%.4f ",exp(-1*PI*PI*(99)*step_t)*sin(PI*step_h*cln1)); }//end of main