int main ( void )
{
setlocale(LC_ALL, "rus");
int I = 10, J = 30, i, j;
double T = 1.0/ pow(3.3, 0.5), h_x = 1.0/ I, h_t = T/ J, epsilon = h_t + pow(h_x, 2), c;
double **u = new double *[I + 1];
for (i = 0; i <= I; i++) u[i] = new double [J + 1];
cout<< "Схема может быть неустойчива при значениях Х :n";
for (i = 0; i <= I; i++)
{
c = 3 * (1.1 - 0.5 * h_x * i) * h_t * pow(h_x, -2);
if (c < 0.5) cout << i * h_x << " ";
}
cout <<"n";
//нулевой слой (j = 0)
for (i = 0; i <= I; i++)
{
u [i][0] = 0.01 * (1 - i * h_x) * i * h_x;
//u [i][0] = 1 - i * h_x; //НУ, несоответствие ГУ и НУ!
}
//последующие слои
for (j = 0; j <= J; j++)
{
for (i = 1; i < I; i++) //расчёт j + 1 - го слоя по j-му
{
u [0][j + 1] = 0; //ГУ u [0][j + 1] = 1;
u [I][j + 1] = 0; //ГУ
u [i][j + 1] = u [i][j] + h_t * (3 * (1.1 - 0.5 * h_x * i) * (u [i + 1][j] -2 * u [i][j] + u [i - 1][j])/ pow(h_x, 2) + exp(h_t * j) - 1);
}
}
int Jv = J/10;
ofstream out;
out.open ("D:\proga7.txt");
out << "U = U(0.6, t):n";
cout << "U = U(0.6, t):n";
for (i = 0; i <= J; i++)
{
out << h_t * i <<"t"<< u [6][i] <<"n";
cout << h_t * i <<"t"<< u [6][i] <<"n";
}
out << "n U = U(x, 0.33):n";
cout << "n U = U(x, 0.33):n";
for (i = 0; i <= I; i++)
{
out << h_x * i <<"t"<< u [i][Jv] <<"n";
cout << h_x * i <<"t"<< u [i][Jv] <<"n";
}
out << "n U = U(x, 0.66):n";
cout << "n U = U(x, 0.66):n";
for (i = 0; i <= I; i++)
{
out << h_x * i <<"t"<< u [i][Jv * 2] <<"n";
cout << h_x * i <<"t"<< u [i][Jv * 2] <<"n";
}
out << "n U = U(x, 1.32):n";
cout << "n U = U(x, 1.32):n";
for (i = 0; i <= I; i++)
{
out << h_x * i <<"t"<< u [i][Jv * 4] <<"n";
cout << h_x * i <<"t"<< u [i][Jv * 4] <<"n";
}
out.close();
getch();
return 0;
}
Метод Конечных Разностей для уравнения в частных производных
Здравствуйте, требуется ваша помощь в нахождении ошибки в программе. Вычисляю по явной схеме. Выдает какие-то очень большие числа.
Задача такая: (формулы у вас не могу ввести корректно.)
[math]U_t=3(1,1-0,5x)U_{xx}+e^t-1[/math]
[math]$ U(0,t)=0$[/math]
[math]$ U(1,t)=0$[/math]
[math]$U(x,0)=0.01(1-x)x$[/math]
Решение нужно найти с точностью [math]$0.0001$[/math] на отрезке [math]$T=1/a^*, где a^*=max a(x,t)$[/math]
Построить графики функций [math]$u(x^*,t), u(x,jt^*)$[/math] где [math]$x^*=0.6, t^*=T/10, j=1,2,4$[/math]
явная разностная схема такая:
([math]$u_t^{j+1}-u_i^j)/tau=3(1,1-0,5x_i)(u_{i+1}^{j}-2u_i^j+u_{i-1}^j)/h^2+e^{t_j}+1$$[/math]
код программы:
выводит следующее:
Добавил также выгрузку динамического массива:
for (i=0; i<I;i++)
{
delete [] u;
}
delete [] u;
ругается:
Задача такая: (формулы у вас не могу ввести корректно.)
[math]U_t=3(1,1-0,5x)U_{xx}+e^t-1[/math]
[math]$ U(0,t)=0$[/math]
[math]$ U(1,t)=0$[/math]
[math]$U(x,0)=0.01(1-x)x$[/math]
Решение нужно найти с точностью [math]$0.0001$[/math] на отрезке [math]$T=1/a^*, где a^*=max a(x,t)$[/math]
Построить графики функций [math]$u(x^*,t), u(x,jt^*)$[/math] где [math]$x^*=0.6, t^*=T/10, j=1,2,4$[/math]
явная разностная схема такая:
([math]$u_t^{j+1}-u_i^j)/tau=3(1,1-0,5x_i)(u_{i+1}^{j}-2u_i^j+u_{i-1}^j)/h^2+e^{t_j}+1$$[/math]
код программы:
Код:
Добавил также выгрузку динамического массива:
for (i=0; i<I;i++)
{
delete [] u;
}
delete [] u;
ругается:
ОС Windows инициировала точку останова в proga7Метод конечных разностей.exe.
Это может быть вызвано повреждением кучи и указывает на ошибку в proga7Метод конечных разностей.exe или в одной из загруженных им DLL.
Возможной причиной так же может быть нажатие пользователем клавиши F12, когда фокус принадлежит proga7Метод конечных разностей.exe
Выведенное на экран окно содержит дополнительные данные для диагностики ошибки
Это может быть вызвано повреждением кучи и указывает на ошибку в proga7Метод конечных разностей.exe или в одной из загруженных им DLL.
Возможной причиной так же может быть нажатие пользователем клавиши F12, когда фокус принадлежит proga7Метод конечных разностей.exe
Выведенное на экран окно содержит дополнительные данные для диагностики ошибки
Код:
#include <iostream>
#include <fstream>
#include <clocale>
#include <math.h>
#include "stdlib.h"
using namespace std;
//------------------------------------------------------------------
/*
dU 1 dU t
---- = 3(1.1 - ---x) ------- + e + 1
dt 2 dx*dx
U(0,t) = 0
U(1,t) = 0
U(x,0) = 0.01(1-x)x = 0.01*x - 0.01*x*x
*/
//------------------------------------------------------------------
//Функция возвращающая коэффициент теплопроводности (сигма)
double Sigma(const double &t,const double &x){
return (3.0 * (1.1 - 0.5 * x));
}
//------------------------------------------------------------------
//Возмущающий член в правой части ДУЧПа
double Funcs(const double &t,const double &x){
return (exp(t) - 1.0);
}
//------------------------------------------------------------------
//Функция возвращающая левое граничное условие
double LeftBC(const double &t,const double &x){
return 0.0;
}
//------------------------------------------------------------------
//Функция возвращающая правое граничное условие
double RightBC(const double &t,const double &x){
return 0.0;
}
//------------------------------------------------------------------
int main(int argc, char *argv[]) {
setlocale(LC_ALL,"Rus");
double length = 1.0, //Протяженность пространства
time, //Время расчета
precision = 1.0,//0.0001,//Точность расчета
dx, //Шаг по пространству
dt, //Шаг по времени
**U, //Расчетный массив
A; //Отношение dt/(dx*dx)
int nX = 10,//Число точек разбиения по пространству
nT = 30;//Число точек разбиения по времени
//Расчет шага по пространству
dx = length/(double)(nX - 1);
//Расчет шага по времени для достижения устойчивости схемы
double Max = 0.0, s;
for(int x=0;x<nX;x++){
s = Sigma(0,dx*x);
Max = (Max > s) ? Max : s;
}
dt = dx*dx*precision/(2.0*Max);
A = dt/(dx*dx);
time = 20.0*dt;//1.0/pow(3.3,0.5);//0.000004;//1.0/Max;
nT = (int)(time/dt + 1.0);
//Выделение памяти под расчетный массив U[time][x]
U = new double *[nT];
for(int t=0;t<nT;t++)
U[t] = new double [nX];
//Инициализация начальных условий
for(int x=0;x<nX;x++)
U[0][x] = 0.01*(1.0 - dx*x)*dx*x;
//Граничные условия для стартового временного слоя
U[0][0] = LeftBC(0,0);
U[0][nX-1] = RightBC(0,dt*(nX-1));
//Вспомогательные переменные
double S, F, T, X;
//Цикл расчета по времени
for(int t=0;t<nT-1;t++){
//Цикл расчета по пространству
for(int x=0;x<nX;x++){
T = dt*t; X = dx*x;
S = Sigma(T,X); F = Funcs(T,X);
//Если на левой границе
if(x == 0){
U[t+1][x] = LeftBC(T,X);
//A*S*(LeftBC(T,X) + U[t][x+1]) + (1.0 - 2.0*A*S)*U[t][x] + dt*F;
} //else//Если на правой границе
if(x == nX-1){
U[t+1][x] = RightBC(T,X);
//A*S*(U[t][x-1] + RightBC(T,X)) + (1.0 - 2.0*A*S)*U[t][x] + dt*F;
}// else {//Если в промежуточных точках
if(x > 0 && x < nX-1){
U[t+1][x] = A*S*(U[t][x-1] + U[t][x+1]) + (1.0 - 2.0*A*S)*U[t][x] + dt*F;
}
}//Конец цикла расчета по пространству
}//Конец цикла расчета по времени
//-------------------------
//Вывод результатов
ofstream out;
out.open ("proga7.txt");
out.precision(10);
out << fixed;
cout.precision(10);
cout << fixed;
//-------------------------
cout << "nX = " << nX << endl;
cout << "nT = " << nT << endl;
cout << "dt = " << dt << endl;
cout << "time = " << time << endl;
cout << "A = " << A << endl << endl;
//-------------------------
out << "nX = " << nX << endl;
out << "nT = " << nT << endl;
out << "dt = " << dt << endl;
out << "time = " << time << endl;
out << "A = " << A << endl << endl;
//-------------------------
int xx = (int)(0.6*(nX-1)/length);
cout << "U = U(0.6,t)" << endl;
out << "U = U(0.6,t)" << endl;
for(int t=0;t<nT;t++){
cout << dt*t << 't' << U[t][xx] << endl;
out << dt*t << 't' << U[t][xx] << endl;
}
cout << endl;
out << endl;
//-------------------------
double tt = time/10.0;
//-------------------------
xx = (int)(tt*(nT-1)*1/time);
cout << "U = U(x," << tt*1 << ")" << endl;
out << "U = U(x," << tt*1 << ")" << endl;
for(int x=0;x<nX;x++){
cout << dx*x << 't' << U[xx][x] << endl;
out << dx*x << 't' << U[xx][x] << endl;
}
cout << endl;
out << endl;
//-------------------------
xx = (int)(tt*(nT-1)*2/time);
cout << "U = U(x," << tt*2 << ")" << endl;
out << "U = U(x," << tt*2 << ")" << endl;
for(int x=0;x<nX;x++){
cout << dx*x << 't' << U[xx][x] << endl;
out << dx*x << 't' << U[xx][x] << endl;
}
cout << endl;
out << endl;
//-------------------------
xx = (int)(tt*(nT-1)*4/time);
cout << "U = U(x," << tt*4 << ")" << endl;
out << "U = U(x," << tt*4 << ")" << endl;
for(int x=0;x<nX;x++){
cout << dx*x << 't' << U[xx][x] << endl;
out << dx*x << 't' << U[xx][x] << endl;
}
//-------------------------
out.close();
//-------------------------
//Освобождение памяти выделенной под массив U
for(int t=0;t<nT;t++){
delete [] U[t];
U[t] = NULL;
}
delete [] U;
U = NULL;
//Конец программы
//-------------------------
system("pause");
return 0;
}
//------------------------------------------------------------------
#include <fstream>
#include <clocale>
#include <math.h>
#include "stdlib.h"
using namespace std;
//------------------------------------------------------------------
/*
dU 1 dU t
---- = 3(1.1 - ---x) ------- + e + 1
dt 2 dx*dx
U(0,t) = 0
U(1,t) = 0
U(x,0) = 0.01(1-x)x = 0.01*x - 0.01*x*x
*/
//------------------------------------------------------------------
//Функция возвращающая коэффициент теплопроводности (сигма)
double Sigma(const double &t,const double &x){
return (3.0 * (1.1 - 0.5 * x));
}
//------------------------------------------------------------------
//Возмущающий член в правой части ДУЧПа
double Funcs(const double &t,const double &x){
return (exp(t) - 1.0);
}
//------------------------------------------------------------------
//Функция возвращающая левое граничное условие
double LeftBC(const double &t,const double &x){
return 0.0;
}
//------------------------------------------------------------------
//Функция возвращающая правое граничное условие
double RightBC(const double &t,const double &x){
return 0.0;
}
//------------------------------------------------------------------
int main(int argc, char *argv[]) {
setlocale(LC_ALL,"Rus");
double length = 1.0, //Протяженность пространства
time, //Время расчета
precision = 1.0,//0.0001,//Точность расчета
dx, //Шаг по пространству
dt, //Шаг по времени
**U, //Расчетный массив
A; //Отношение dt/(dx*dx)
int nX = 10,//Число точек разбиения по пространству
nT = 30;//Число точек разбиения по времени
//Расчет шага по пространству
dx = length/(double)(nX - 1);
//Расчет шага по времени для достижения устойчивости схемы
double Max = 0.0, s;
for(int x=0;x<nX;x++){
s = Sigma(0,dx*x);
Max = (Max > s) ? Max : s;
}
dt = dx*dx*precision/(2.0*Max);
A = dt/(dx*dx);
time = 20.0*dt;//1.0/pow(3.3,0.5);//0.000004;//1.0/Max;
nT = (int)(time/dt + 1.0);
//Выделение памяти под расчетный массив U[time][x]
U = new double *[nT];
for(int t=0;t<nT;t++)
U[t] = new double [nX];
//Инициализация начальных условий
for(int x=0;x<nX;x++)
U[0][x] = 0.01*(1.0 - dx*x)*dx*x;
//Граничные условия для стартового временного слоя
U[0][0] = LeftBC(0,0);
U[0][nX-1] = RightBC(0,dt*(nX-1));
//Вспомогательные переменные
double S, F, T, X;
//Цикл расчета по времени
for(int t=0;t<nT-1;t++){
//Цикл расчета по пространству
for(int x=0;x<nX;x++){
T = dt*t; X = dx*x;
S = Sigma(T,X); F = Funcs(T,X);
//Если на левой границе
if(x == 0){
U[t+1][x] = LeftBC(T,X);
//A*S*(LeftBC(T,X) + U[t][x+1]) + (1.0 - 2.0*A*S)*U[t][x] + dt*F;
} //else//Если на правой границе
if(x == nX-1){
U[t+1][x] = RightBC(T,X);
//A*S*(U[t][x-1] + RightBC(T,X)) + (1.0 - 2.0*A*S)*U[t][x] + dt*F;
}// else {//Если в промежуточных точках
if(x > 0 && x < nX-1){
U[t+1][x] = A*S*(U[t][x-1] + U[t][x+1]) + (1.0 - 2.0*A*S)*U[t][x] + dt*F;
}
}//Конец цикла расчета по пространству
}//Конец цикла расчета по времени
//-------------------------
//Вывод результатов
ofstream out;
out.open ("proga7.txt");
out.precision(10);
out << fixed;
cout.precision(10);
cout << fixed;
//-------------------------
cout << "nX = " << nX << endl;
cout << "nT = " << nT << endl;
cout << "dt = " << dt << endl;
cout << "time = " << time << endl;
cout << "A = " << A << endl << endl;
//-------------------------
out << "nX = " << nX << endl;
out << "nT = " << nT << endl;
out << "dt = " << dt << endl;
out << "time = " << time << endl;
out << "A = " << A << endl << endl;
//-------------------------
int xx = (int)(0.6*(nX-1)/length);
cout << "U = U(0.6,t)" << endl;
out << "U = U(0.6,t)" << endl;
for(int t=0;t<nT;t++){
cout << dt*t << 't' << U[t][xx] << endl;
out << dt*t << 't' << U[t][xx] << endl;
}
cout << endl;
out << endl;
//-------------------------
double tt = time/10.0;
//-------------------------
xx = (int)(tt*(nT-1)*1/time);
cout << "U = U(x," << tt*1 << ")" << endl;
out << "U = U(x," << tt*1 << ")" << endl;
for(int x=0;x<nX;x++){
cout << dx*x << 't' << U[xx][x] << endl;
out << dx*x << 't' << U[xx][x] << endl;
}
cout << endl;
out << endl;
//-------------------------
xx = (int)(tt*(nT-1)*2/time);
cout << "U = U(x," << tt*2 << ")" << endl;
out << "U = U(x," << tt*2 << ")" << endl;
for(int x=0;x<nX;x++){
cout << dx*x << 't' << U[xx][x] << endl;
out << dx*x << 't' << U[xx][x] << endl;
}
cout << endl;
out << endl;
//-------------------------
xx = (int)(tt*(nT-1)*4/time);
cout << "U = U(x," << tt*4 << ")" << endl;
out << "U = U(x," << tt*4 << ")" << endl;
for(int x=0;x<nX;x++){
cout << dx*x << 't' << U[xx][x] << endl;
out << dx*x << 't' << U[xx][x] << endl;
}
//-------------------------
out.close();
//-------------------------
//Освобождение памяти выделенной под массив U
for(int t=0;t<nT;t++){
delete [] U[t];
U[t] = NULL;
}
delete [] U;
U = NULL;
//Конец программы
//-------------------------
system("pause");
return 0;
}
//------------------------------------------------------------------
Скинь, как нибудь, исходные формулы в читаемой форме (рисунком, хотя бы).
да скинь, а то так не понятно
Meander
Большое вам, человеческое спасибо! чуть не плачу от счастья! уже не надеялся что кто нибудь поможет!!!
Большое вам, человеческое спасибо! чуть не плачу от счастья! уже не надеялся что кто нибудь поможет!!!