Конечно-разностные схемы моделирования распространения волн
Дипломная работа - Компьютеры, программирование
Другие дипломы по предмету Компьютеры, программирование
um; j++){(int k = 0; k < x_num; k++){(fout, "%f\t", out_array[i][j][k]);
}(fout, "\n");
}(fout, "\n\n");
}(fout);
}:void setInitialConditions() = 0;void setBoundaryConditions(int n) = 0;void calculate() = 0;(int N, int M, double Dx, double H_temp){= Dx;= H_temp;= dx/100;_num = N;_num = M;= new double**[t_num];= new double**[t_num];= new double**[t_num];(int i = 0; i < t_num; i++){[i] = new double*[x_num];[i] = new double*[x_num];[i] = new double*[x_num];(int j = 0; j < x_num; j++){[i][j] = new double[x_num];[i][j] = new double[x_num];[i][j] = new double[x_num];
}
}= new double*[x_num];(int i = 0; i < x_num; i++){[i] = new double[x_num];
}
}() {}
~Scheme(){(int i = 0; i < t_num; i++){(int j = 0; j < x_num; j++){[] u[i][j];[] v[i][j];[] eta[i][j];
}[] u[i];[] v[i];[] eta[i];
}[] u;[] v;[] eta;(int i = 0; i < x_num; i++){[] etaTemp[i];
}[] etaTemp;
}clear(){(int i = 0; i < t_num; i++){(int j = 0; j < x_num; j++){(int k = 0; k < x_num; k++){[i][j][k] = 0.0;[i][j][k] = 0.0;[i][j][k] = 0.0;
}
}
}(int j = 0; j < x_num; j++){(int k = 0; k < x_num; k++){[j][k] = 0.0;
}
}
}ouput(){* foutConfig = fopen("config.txt", "w");(foutConfig, "%d\n%lf\n", x_num, dx);(foutConfig, "%d\n%lf\n", t_num, dt);(foutConfig);
("eta.txt", eta);
}
};SchemeN: public Scheme {* a;* b;* c;* f;* u_temp;* v_temp;* d;* e;:void setInitialConditions(){x0 = x_num/2;y0 = x_num/2;eta0 = 1.0;a = 0.2;(int i = 0; i < x_num; i++){(int j = 0; j < x_num; j++){[0][i][j] = 0.0;[0][i][j] = 0.0;
//eta[0][i][j] = H + eta0*exp( - a*a*((dx*i-x0)*(dx*i-x0) + (dx*j-y0)*(dx*j-y0)));[0][i][j]= exp(-dx*(i-x0)*3*dx*(i-x0)-dx*(j-y0)*dx*3*(j-y0));[1][i][j] = eta[0][i][j];
}
}
}void setBoundaryConditions(int n){(int i = 1; i < x_num-1; i++){[n][0][i] = 2*eta[n][1][i] - eta[n][2][i];[n][x_num-1][i] = 2*eta[n][x_num-2][i] - eta[n][x_num-3][i];[n][i][0] = 2*eta[n][i][1] - eta[n][i][2];[n][i][x_num-1] = 2*eta[n][i][x_num-2] - eta[n][i][x_num-3];[n][0][i] = 2*u[n][1][i] - u[n][2][i];[n][x_num-1][i] = 2*u[n][x_num-2][i] - u[n][x_num-3][i];[n][i][0] = 2*u[n][i][1] - u[n][i][2];[n][i][x_num-1] = 2*u[n][i][x_num-2] - u[n][i][x_num-3];[n][0][i] = 2*v[n][1][i] - v[n][2][i];[n][x_num-1][i] = 2*v[n][x_num-2][i] - v[n][x_num-3][i];[n][i][0] = 2*v[n][i][1] - v[n][i][2];[n][i][x_num-1] = 2*v[n][i][x_num-2] - v[n][i][x_num-3];
}[n][0][0] = 2*eta[n][1][0] - eta[n][2][0];[n][x_num-1][0] = 2*eta[n][x_num-1][1] - eta[n][x_num-1][2];[n][0][x_num-1] = 2*eta[n][1][x_num-1] - eta[n][2][x_num-1];[n][x_num-1][x_num-1] = 2*eta[n][x_num-1][x_num-2] - eta[n][x_num-1][x_num-3];[n][0][0] = 2*u[n][1][0] - u[n][2][0];[n][x_num-1][0] = 2*u[n][x_num-1][1] - u[n][x_num-1][2];[n][0][x_num-1] = 2*u[n][1][x_num-1] - u[n][2][x_num-1];[n][x_num-1][x_num-1] = 2*u[n][x_num-1][x_num-2] - u[n][x_num-1][x_num-3];[n][0][0] = 2*v[n][1][0] - v[n][2][0];[n][x_num-1][0] = 2*v[n][x_num-1][1] - v[n][x_num-1][2];[n][0][x_num-1] = 2*v[n][1][x_num-1] - v[n][2][x_num-1];[n][x_num-1][x_num-1] = 2*v[n][x_num-1][x_num-2] - v[n][x_num-1][x_num-3];
}solveMatrix(double * a, double * b, double *c, double * f, double * x){[1] = -b[0]/c[0];[1] = f[0]/c[0];(int i = 1; i < x_num-1; ++i) {[i+1] = -b[i] / (a[i]*d[i]+c[i]);[i+1] = (f[i] - a[i]*e[i]) / (c[i]+a[i]*d[i]);
}[x_num-1] = (f[x_num-1] - a[x_num-1]*e[x_num-1]) / (c[x_num-1]+a[x_num-1]*d[x_num-1]);(int i = x_num-2; i > 0; --i) {[i] = d[i+1]*x[i+1]+e[i+1];
}
}prepareEtaTemp(int n){(int i = 1; i < x_num-1; i++){(int j = 1; j < x_num-1; j++){[i][j] = eta[n][i][j] - dt*((H + eta[n][i][j])*(u[n][i+1][j] - u[n][i-1][j])/(2.0*dx)
+ (H + eta[n][i][j])*(v[n][i][j+1] - u[n][i][j-1])/(2.0*dx)
+ u[n][i][j]*(eta[n][i+1][j] - eta[n][i-1][j])/(2.0*dx));
}
}
}calculateEta(int n){(int i = 1; i < x_num-1; i++){(int j = 1; j < x_num-1; j++){[n+1][i][j] = eta[n][i][j] - dt*((H + eta[n][i][j])*(u[n+1][i+1][j] - u[n+1][i-1][j] + u[n][i+1][j] - u[n][i-1][j])/dx
+ (H + eta[n][i][j])*(v[n+1][i][j+1] - v[n+1][i][j-1] + v[n][i][j+1] - v[n][i][j-1])/dx
+ ((v[n+1][i][j] + v[n][i][j])/2.0)*((eta[n][i][j] - eta[n][i][j-1])/(2*dx)));
}
}
}calculateVx(int n){(int j = 1; j < x_num-1; j++){[0] = a[x_num-1] = 0;[0] = b[x_num-1] = 0;[0] = c[x_num-1] = 1;(int i = 0; i < x_num; i++) {_temp[i] = v[n][i][j];
}[0] = v_temp[0];[x_num-1] = v_temp[x_num-1];(int i = 1; i < x_num-1; i++) {[i] = -u[n][i][j]/(4.0*dx);[i] = 1.0/dt;[i] = u[n][i][j]/(4.0*dx);[i] = -v[n][i][j]/(4.0*dx)*v[n+1][i][j+1] + v[n][i][j]/(4.0*dx)*v[n+1][i][j-1]
+ v[n][i][j]/dt-u[n][i][j]*(v[n][i+1][j] - v[n][i-1][j])/(4.0*dx)
v[n][i][j]*(v[n][i][j+1] - v[n][i][j-1])/(4.0*dx)
g*(etaTemp[i][j+1] - etaTemp[i][j-1] + eta[n][i][j+1] - eta[n][i][j-1])/(4.0*dx);
}(a, b, c, f, v_temp);(int i = 0; i < x_num; i++) {[n+1][i][j] = v_temp[i];
}
}
}calculateVy(int n){(int i = 1; i < x_num-1; i++){[0] = a[x_num-1] = 0;[0] = b[x_num-1] = 0;[0] = c[x_num-1] = 1;(int j = 0; j < x_num; j++) {_temp[j] = v[n][i][j];
}[0] = v_temp[0];[x_num-1] = v_temp[x_num-1];(int j = 1; j < x_num-1; j++) {[j] = -v[n][i][j]/(4.0*dx);[j] = 1.0/dt;[j] = v[n][i][j]/(4.0*dx);[j] = (-v[n][i][j]/(4.0*dx))*v[n+1][i+1][j] + (u[n][i][j]/(4.0*dx))*v[n+1][i-1][j]
+ v[n][i][j]/dt - u[n][i][j]*(v[n][i+1][j]-v[n][i-1][j])/(4.0*dx)
v[n][i][j]*(v[n][i][j+1] - v[n][i][j-1])/(4.0*dx)
g*(etaTemp[i][j+1] - etaTemp[i][j-1] + eta[n][i][j+1] - eta[n][i][j-1])/(4.0*dx);
}(a, b, c, f, v_temp);(int j = 0; j < x_num; j++) {[n+1][i][j] = v_temp[j];
}
}
}calculateUx(int n){(int j = 1; j < x_num-1; j++){[0] = a[x_num-1] = 0;[0] = b[x_num-1] = 0;[0] = c[x_num-1] = 1;(int i = 0; i < x_num; i++) {_temp[i] = u[n][i][j];
}[0] = u_temp[0];[x_num-1] = u_temp[x_num-1];(int i = 1; i < x_num-1; i++) {[i] = -u[n][i][j]/(4.0*dx);[i] = 1.0/dt;[i] = u[n][i][j]/(4.0*dx);[i] = (-v[n][i][j]/(4.0*dx))*u[n+1][i][j+1] + (v[n][i][j]/(4.0*dx))*u[n+1][i][j-1]
+ u[n][i][j]/dt - u[n][i][j]*(u[n][i+1][j] - u[n][i-1][j])/(4.0*dx)
v[n][i][j]*(u[n][i][j+1] - u[n][i][j-1])/(4.0*dx)
g*(etaTemp[i+1][j] - etaTemp[i-1][j] + eta[n][i+1][j] - eta[n][i-1][j])/(4.0*dx);
}(a, b, c, f, u_temp);(int i = 0; i < x_num; i++) {[n+1][i][j] = u_temp[i];
}
}
}calculateUy(int n){(int i = 1; i < x_num-1; i++){[0] = a[x_num-1] = 0;[0] = b[x_num-1] = 0;[0] = c[x_num-1] = 1;(int j = 0; j < x_num; j++) {_temp[j] = u[n][i][j];
}[0] = u_temp[0];[x_num-1] = u_temp[x_num-1];(int j = 1; j < x_num-1; j++) {[j] = -v[n][i][j]/(4.0*dx);[j] = 1.0/dt;[j] = v[n][i][j]/(4.0*dx);[j] = (-u[n][i][j]/(4.0*dx))*u[n+1][i+1][j] + (u[n][i][j]/(4.0*dx))*u[n+1][i-1][j]
+ u[n][i][j]/dt - u[n][i][j]*(u[n][i+1][j] - u[n][i-1][j])/(4.0*dx)
v[n][i][j]*(u[n][i][j+1] - u[n][i][j-1])/(4.0*dx)
g*(etaTemp[i+1][j] - etaTemp[i-1][j] + eta[n][i+1][j] - eta[n][i-1][j])/(4.0*dx);
}(a, b, c, f, u_temp);(int j = 0; j < x_num; j++) {[n+1][i][j] = u_temp[j];
}
}
}void calculate(){(int i = 0; i < t_num-1; i++){(i);(i);(i);(i);(i);(i);(i);
}
}(int N, int M, double Dx, double H_temp){= Dx;= H_temp;= dx/100;_num = N;_num = M;= new double**[t_num];= new double**[t_num];= new double**[t_num];(int i = 0; i < t_num; i++){[i] = new double*[x_num];[i] = new double*[x_num];[i] = new double*[x_num];(int j = 0; j < x_num; j++){[i][j] = new double[x_num];[i][j] = new double[x_num];[i][j] = new double[x_num];
}
}
= new double*[x_num];(int i = 0; i < x_num; i++){[i] = new double[x_num];
}= new double[x_num];= new double[x_num];= new double[x_num];= new double[x_num];= new double[x_num];= new double[x_num];_temp = new double[x_num];_temp = new double[x_num];
}
};main(){test(200, 50, 0.1, 1.0);.clear();.setInitialConditions();.calculate();.ouput();("pause");1;
}
3.7.3Исходный код программы для вывода в среде MATLAB
Представлен исходный код программы для вывода в среде MATLAB.
код программа дифференциальный уравнение
Исходный код программы для вывода в среде MATLAB
clear;config.txt;= config(2);= config(4);= 0:dx:(config(1)-1)*dx;= x;= config(3);= moviein(num);= 1;eta.txt;= 0;= dx*(config(1)-1);r = 1:num= eta(n:n+config(1)-1,:);(x,y,z);([0 xmax 0 xmax -0 2]);('x');('y');('z');(:,r) = getframe;= n+config(1);= 0;= 10;(M, repeat, fps);
3.8Тестовые примеры для программы, реализующей явную схему
Программа реализованная на C/C++ в качестве выходных данных имеет 2 файла: config.txt - содержащий информацию о расчётной сетке и u.txt - содержащую распределение температур. Примеры выходных данных можно посмотреть в прикреплённых файлах.
На рисунках 12 - 13 изображены результаты работы программы строящей графики в среде MATLAB в разные моменты времени.
Рисунок 12 - состояние волны в начальный момент времени
Рисунок 13 - с