Crank_Nicolson.h
#include <GL/glut.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <malloc.h>
#include <math.h>
int xnseg_ = 100;
int ynseg_ = 100;
int tnseg_ = 100;
struct Heat {
double u[101][101];
};
struct Heat heat[101];
struct Heat_h {
double u_h[101][101];
};
struct Heat_h heat_h[100];
double xx[101];
double yy[101];
/////////////////////////////////////////////////////////////////////////
class Crank_Nicolson {
public :
Crank_Nicolson(void);
void clear(void);
void clear_index(void);
~Crank_Nicolson(void) { clear(); }
public :
void set_init_funct0(double (*_funct0)(double _x, double _y)) { funct0_ = _funct0; }
void set_init_funcx0(double (*_funcx0)(double _y, double _t)) { funcx0_ = _funcx0; } //初始条件 u(x,0)=f(x)
void set_init_funcxn(double (*_funcxn)(double _y, double _t)) { funcxn_ = _funcxn; }
void set_init_funcy0(double (*_funcy0)(double _x, double _t)) { funcy0_ = _funcy0; } //初始条件 u(x,0)=f(x)
void set_init_funcyn(double (*_funcyn)(double _x, double _t)) { funcyn_ = _funcyn; }
void setSegment(double _xbgn, double _xend, double _ybgn, double _yend); //_xbgn,_xend分别为数组 x 的初值和末值; _xnseg为将x 分成的段数
void init_arrays(double _tbgn, double _tvary); //_tbgn为起始时间, _tseg为将时间分成的段数, _tvary为相邻时间的差值
void set_c(double _c); //设置温度系数c
void get_rx(void); //计算系数r
void get_ry(void);
void init_index_hf(int layer, int row);
void guassi_hf();
void calculate_hf_onerow(int layer, int row);
void calculate_hf_onelayer(int layer);
void init_index(int layer, int col);
void guassi();
void calculate_onecol(int layer, int col);
void calculate_onelayer(int layer);
void calculate(void); //计算不同时间不同位置处的温度, 并将结果放入数组二维u
void show_layer(int layer);
private :
double xbgn_, xend_, ybgn_, yend_, tbgn_, tvary_; //tvary_为相邻时间的差值
double c_, rx_, ry_;
double (*funct0_)(double,double);
double (*funcx0_)(double,double);
double (*funcxn_)(double,double);
double (*funcy0_)(double,double);
double (*funcyn_)(double,double);
double *y, *x, *t; //t 存放时间, x存放位置
double *t_h;
double *index1;
double *index2;
double *index3;
double *result;
};
file.h
/////////////////////////////////////////////////////////////////////////////////////////////
/////////////////////////////////////////////////////////////////////////////////////////////
#include "Crank_Nicolson.h"
/////////////////////////////////////////////////////////////////////////////////////////
Crank_Nicolson::Crank_Nicolson() {
xbgn_ = 0; xend_ = 0; ybgn_ = 0; yend_ = 0;
tbgn_ = 0; tvary_ = 0;
c_ = 0; rx_ = 0; ry_ = 0;
y = NULL; x = NULL; t = NULL; t_h = NULL;
}
////////////////////////////////////////////////////////////////////////////////////////////
void Crank_Nicolson::clear() {
if(y)
delete []y;
if(x)
delete []x;
if(t)
delete []t;
}
/////////////////////////////////////////////////////////////////////////////////
void Crank_Nicolson::setSegment(double _xbgn, double _xend, double _ybgn, double _yend) {
if(_xbgn >= _xend) {
std::cout << "Err: range of X is wrong." << '\n';
exit(0);
}
if(_ybgn >= _yend) {
std::cout << "Err: range of y is wrong." << '\n';
exit(0);
}
xbgn_ = _xbgn;
xend_ = _xend;
ybgn_ = _ybgn;
yend_ = _yend;
}
//////////////////////////////////////////////////////////////////////////////
void Crank_Nicolson::init_arrays(double _tbgn, double _tvary) {
clear();
if(_tbgn < 0) {
std::cout << "Err: t is NOT defined at all." << '\n';
exit(0);
}
if(_tvary < 0) {
std::cout << "Err: tvary is NOT defined at all." << '\n';
exit(0);
}
tbgn_ = _tbgn; tvary_ = _tvary;
index1 = (double *) malloc ((xnseg_-2)*sizeof(double));
index2 = (double *) malloc ((xnseg_-1)*sizeof(double));
index3 = (double *) malloc ((xnseg_-2)*sizeof(double));
result = (double *) malloc ((xnseg_-1)*sizeof(double));
y = (double *)malloc((ynseg_+1)*sizeof(double));
x = (double *)malloc((xnseg_+1)*sizeof(double));
t = (double *)malloc((tnseg_+1)*sizeof(double));
t_h = (double *)malloc((tnseg_)*sizeof(double));
*(t+0) = tbgn_;
for(int i = 1; i <= tnseg_; i++)
*(t+i) = *(t+i-1) + tvary_;
*(t_h+0) = tbgn_+(1.0/2)*tvary_;
for(i = 1; i < tnseg_; i++)
*(t_h+i) = *(t_h+i-1) + tvary_;
*(y+0) = ybgn_;
yy[0]=*(y+0);
double yi = (yend_ - ybgn_) / ynseg_;
for(i = 1; i <= ynseg_; i++) {
*(y+i) = *(y+i-1) + yi;
yy[i]=*(y+i);
}
*(x+0) = xbgn_;
xx[0]=*(x+0);
double xi = (xend_ - xbgn_) / xnseg_;
for (i = 1; i <= xnseg_; i++) {
*(x+i) = *(x+ i-1) + xi;
xx[i]=*(x+i);
}
for (int j=0; j<=tnseg_; j++) {
for (i = 0; i <= ynseg_ ; i++) {
heat[j].u[0][i] = funcx0_(*(y+i),*(t+i));
heat[j].u[xnseg_][i] = funcxn_(*(y+i),*(t+i));
}
for (i = 1; i < xnseg_ ; i++) {
heat[j].u[i][0] = funcy0_(*(x+i),*(t+i));
heat[j].u[i][ynseg_] = funcyn_(*(x+i),*(t+i));
}
}
for (j=0; j<tnseg_; j++) {
for (i = 0; i <= ynseg_ ; i++) {
heat_h[j].u_h[0][i] = funcx0_(*(y+i),*(t_h+i));
heat_h[j].u_h[xnseg_][i] = funcxn_(*(y+i),*(t_h+i));
}
for (i = 1; i < xnseg_ ; i++) {
heat_h[j].u_h[i][0] = funcy0_(*(x+i),*(t_h+i));
heat_h[j].u_h[i][ynseg_] = funcyn_(*(x+i),*(t_h+i));
}
}
//initial condition: u(x,0)=f(x)
for (i = 1; i< xnseg_; i++)
for(int j=1; j<ynseg_; j++)
heat[0].u[i][j] = funct0_ (*(x+i),*(y+j)); //equation func_(float) has been identified before
}
//////////////////////////////////////////
void Crank_Nicolson::set_c(double _c) {
c_ = _c;
}
////////////////////////////////////////////////
void Crank_Nicolson::get_rx() {
double xi = (xend_ - xbgn_) / xnseg_;
rx_ = (c_*tvary_) / (xi*xi);
std::cout << "rx = " << rx_ << '\n';
}
void Crank_Nicolson::get_ry() {
double yi = (yend_ - ybgn_) / ynseg_;
ry_ = (c_*tvary_) / (yi*yi);
std::cout << "ry = " << ry_ << '\n';
}
////////////////////////////////////////////////////////////////////
////////////////////////////////////////new!!!!!!!!
void Crank_Nicolson::init_index_hf(int layer, int row) { //从0 1 开始 row 到 ynseg_-1 结束
*(index2+1) = 1+rx_;
*(index3+1) = -1*rx_*(1.0/2);
*(result+1) = 1.0/2*rx_*heat[layer].u[1][row-1] + (1-rx_)*heat[layer].u[1][row] + 1.0/2*rx_*heat[layer].u[1][row+1] +1.0/2*rx_*heat_h[layer].u_h[0][row];
for(int i = 2; i < xnseg_-1; i++) {
*(index1+i) = -1*rx_*(1.0/2);
*(index2+i) = 1+rx_;
*(index3+i) = -1*rx_*(1.0/2);
*(result+i) = 1.0/2*rx_*heat[layer].u[i][row-1] + (1-rx_)*heat[layer].u[i][row] + 1.0/2*rx_*heat[layer].u[i][row+1];
}
*(index1+i) = -1*rx_*(1.0/2);
*(index2+i) = 1+rx_;
*(result+i) = 1.0/2*rx_*heat[layer].u[i][row-1] + (1-rx_)*heat[layer].u[i][row] + 1.0/2*rx_*heat[layer].u[i][row+1] + 1.0/2*rx_*heat_h[layer].u_h[xnseg_][row];
}
void Crank_Nicolson::guassi_hf() {
for (int i = 2; i < xnseg_-1; i++) {
*(index2+i) = (*(index2+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) ) + (*(index3+i-1));
*(index3+i) = (*(index3+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) );
*(result+i) = (*(result+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) ) + (*(result+i-1));
*(index1+i) = 0;
}
*(index2+i) = (*(index2+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) ) + (*(index3+i-1));
*(result+i) = (*(result+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) ) + (*(result+i-1));
*(index1+i) = 0;
}
void Crank_Nicolson::calculate_hf_onerow(int layer, int row) { //0,1
init_index_hf(layer, row);
guassi_hf();
double y_pre;
y_pre = ( (*(result+xnseg_-1)) ) / (*(index2+xnseg_-1));
double y_next = y_pre;
heat_h[layer].u_h[xnseg_-1][row] = y_next;
for (int i = xnseg_-2; i > 1; i--) {
y_pre = ( (*(result+i)) - ( (*(index3+i)) * y_next ) ) / (*(index2+i));
heat_h[layer].u_h[i][row] = y_pre;
y_next = y_pre;
}
y_pre = ( *(result+1) - (*(index3+1))*y_next) / (*(index2+1));
heat_h[layer].u_h[i][row] = y_pre;
}
void Crank_Nicolson::calculate_hf_onelayer(int layer) { //0~tnseg_-1
int row;
for(row=1; row<ynseg_; row++)
calculate_hf_onerow(layer, row);
}
void Crank_Nicolson::init_index(int layer, int col) { //从0 1 开始 col 到 xnseg_-1 结束
*(index2+1) = 1+rx_;
*(index3+1) = -1*rx_*(1.0/2);
*(result+1) = 1.0/2*rx_*heat_h[layer].u_h[col-1][1] + (1-rx_)*heat_h[layer].u_h[col][1] + 1.0/2*rx_*heat_h[layer].u_h[col+1][1] +(1.0/2)*rx_*heat[layer+1].u[col][0];
for(int i = 2; i < ynseg_-1; i++) {
*(index1+i) = -1*rx_*(1.0/2);
*(index2+i) = 1+rx_;
*(index3+i) = -1*rx_*(1.0/2);
*(result+i) = (1.0/2)*rx_*heat_h[layer].u_h[col-1][i] + (1-rx_)*heat_h[layer].u_h[col][i] + (1.0/2)*rx_*heat_h[layer].u_h[col+1][i];
}
*(index1+i) = -1*rx_*(1.0/2);
*(index2+i) = 1+rx_;
*(result+i) = 1.0/2*rx_*heat_h[layer].u_h[col-1][i] + (1-rx_)*heat_h[layer].u_h[col][i] + 1.0/2*rx_*heat_h[layer].u_h[col+1][i] + 1.0/2*rx_*heat[layer+1].u[col][ynseg_];
}
void Crank_Nicolson::guassi() {
for (int i = 2; i < ynseg_-1; i++) {
*(index2+i) = (*(index2+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) ) + (*(index3+i-1));
*(index3+i) = (*(index3+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) );
*(result+i) = (*(result+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) ) + (*(result+i-1));
*(index1+i) = 0;
}
*(index2+i) = (*(index2+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) ) + (*(index3+i-1));
*(result+i) = (*(result+i)) * ( (-1)*(*(index2+i-1)) / (*(index1+i)) ) + (*(result+i-1));
*(index1+i) = 0;
}
void Crank_Nicolson::calculate_onecol(int layer, int col) { //0,1
init_index(layer, col);
guassi();
double y_pre;
y_pre = ( (*(result+ynseg_-1)) ) / (*(index2+ynseg_-1));
double y_next = y_pre;
heat[layer+1].u[col][ynseg_-1] = y_next;
for (int i = ynseg_-2; i > 1; i--) {
y_pre = ( (*(result+i)) - ( (*(index3+i)) * y_next ) ) / (*(index2+i));
heat[layer+1].u[col][i] = y_pre;
y_next = y_pre;
}
y_pre = ( *(result+1) - (*(index3+1))*y_next) / (*(index2+1));
heat[layer+1].u[col][i] = y_pre;
}
void Crank_Nicolson::calculate_onelayer(int layer) { //0~tnseg_-1
int col;
for(col=1; col<xnseg_; col++)
calculate_onecol(layer, col);
}
///new!!!!!!!!!!!!!!!!!!!!!!!
//////////////////////////////////////////////////////////////////////
void Crank_Nicolson::calculate() {
int layer;
for(layer=0; layer<tnseg_; layer++) {
calculate_hf_onelayer(layer);
calculate_onelayer(layer);
}
}
void Crank_Nicolson::show_layer(int layer) {
std::cout << '\n' << "layer= " << layer << '\n';
for (int i=0; i<=xnseg_; i=i+5) {
std::cout<< '\n' << "x=" << *(x+i) << '\n';
for(int j=0; j<=ynseg_; j=j+10) {
if(heat[layer].u[i][j]<0)
heat[layer].u[i][j]=0;
std::cout<< heat[layer].u[i][j] << " ";
}
}
std::cout << '\n';
}
Crank_Nicolson.cpp
#include <GL/glut.h>
#include <stdio.h>
#include <stdlib.h>
#define PI 3.1415926
void init(void);
void reshape(int w,int h);
void display(void);
void keyboard(unsigned char key, int x, int y);
int layer = 0;
double maxt0 = 100.0;
#include "file.h"
#include <math.h>
Crank_Nicolson rank_Nicolson;
///////////////////////////////////////////
double funct0(double x, double y) {
return 100-sqrt((x-50.0)*(x-50.0)+(y-50.0)*(y-50.0));
// return 100*sin(PI/10*x)+100*sin(PI/100*y);
}
double funcx0(double y, double t) {
return 0.0;
}
double funcxn(double y, double t) {
return 0.0;
}
double funcy0(double x, double t) {
return 0.0;
}
double funcyn(double x, double t) {
return 0.0;
}
//////////////////////////////////////
int main (int argc, char **argv) {
rank_Nicolson.set_init_funct0(funct0);
rank_Nicolson.set_init_funcx0(funcx0);
rank_Nicolson.set_init_funcxn(funcxn);
rank_Nicolson.set_init_funcy0(funcy0);
rank_Nicolson.set_init_funcyn(funcyn);
rank_Nicolson.setSegment(0.0, 100.0, 0.0, 100.0); //步长为1
rank_Nicolson.init_arrays(0.0, 1);
rank_Nicolson.set_c(10);
rank_Nicolson.get_rx();
rank_Nicolson.get_ry();
rank_Nicolson.calculate();
glutInit(&argc, argv);
glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB); //用GLUT_DOUBLE模式可以消除动画程序的闪烁现象
glutInitWindowSize (500, 500);
glutInitWindowPosition (100, 100);
glutCreateWindow("myheat");
init();
glutDisplayFunc(display);
glutReshapeFunc(reshape);
glutKeyboardFunc(keyboard);
glutMainLoop();
return 0;
}
void init (void)
{
glClearColor (1.0, 1.0, 1.0, 0.0);
glShadeModel(GL_SMOOTH);
}
void reshape(int w, int h)
{
glViewport (0, 0, (GLsizei) w, (GLsizei) h); //定义适口大小,(x,y)为左下角坐标,(width,height)为高和宽
glMatrixMode(GL_PROJECTION);
glLoadIdentity();
//gluLookAt(50.0,30.0,50.0,0.0,0.0,-1.0,1.0,0.0,0.0);
glOrtho(-200, 300, 200, -400, -500, 500); //(left,right,bottom,top,near,far)
//创建正交平行视景体矩阵,并把它与当前矩阵相乘,(left,bottom,-near)和(right,top,-near)是近侧裁剪平面上的点,分别映射到适口窗口的左下角和右上角;
//(left,bottom,-far)和(right,top,-far)是远侧裁剪平面上的点,分别映射到视口窗口的左下角和右上角.
glRotatef(90,0,1,0);
glRotatef(90,1,0,0);
glTranslated(10.0,-30,0.0);
glScalef(2.0,2.0,2.0);
glMatrixMode(GL_MODELVIEW);
glLoadIdentity();
}
void display() {
glClear (GL_COLOR_BUFFER_BIT);
GLint i, j;
GLdouble x_lf;
for(i=0; i<xnseg_; i++) {
for(j=0; j<ynseg_; j++) {
glBegin(GL_POINTS);
//glPolygonMode(GL_FRONT_AND_BACK,GL_LINES);
x_lf = heat[layer].u[i][j] / maxt0;
if(x_lf>1) x_lf=1;
glColor3f(0.0+x_lf,0.0,1.0-(0.0+x_lf));
glVertex3f(xx[i],yy[j],heat[layer].u[i][j]);
glEnd();
}
}
glutSwapBuffers(); //与双缓冲模式相对应,若是单缓冲模式,则用glFlush()
rank_Nicolson.show_layer(layer);
}
void keyboard(unsigned char key, int x, int y) {
if(key == 'g' || key == 'G') {
layer++;
display();
}
}