Язык С++ как инструмент моделирования на основе решения дифференциальных уравнений в частных производных - page 9

7
}
//при t начальном U=sin(X):
for(j=0; j<nx; j++){
U[0][j]=sin(X[j]);
}
//расчет в цикле всей поверхности:
int t;
for(t=0; t<nt-1; t++){
for(i=1; i<nx-1; i++){
U[t+1][i]=U[t][i]+dt/dx/dx*a*
(U[t][i+1]-2*U[t][i]+U[t][i-1]);//////+++
}
U[t+1][0]=U[t+1][1]-dx*exp(-a*T[t+1]);////////+++
U[t+1][(int)nx-1]=U[t+1][(int)nx-2]-
dx*exp(-a*T[t+1]);
}
//аналитическое решение:
double** U_an = new double*[(int)nt];
for ( i = 0; i < nt; i++) {
U_an[i] = new double[(int)nx];
}
for(t=0; t<nt; t++){
for(i=0; i<nx; i++){
U_an[t][i]=exp(-a*T[t])*sin(X[i]);
}
}
//среднеквадратическая погрешность:
double error,er1=0,er2=0;
for(t=0; t<nt; t++){
for(i=0; i<nx; i++){
er1=er1+pow(U[t][i]-U_an[t][i],2.0);
er2=er2+pow(U[t][i],2.0);
}
}
error=pow(er1/er2,0.5);
//вывод в файл MatLAB:
fstream m;
m.open("shema1 10 10.m", ios::out);
m<<"% задача 1 схема 1 размер 10 10"<<endl;
m<<"T = [ ";
for(t=0; t<nt; t++){
m<<T[t]<<" ";
}
m<<"];"<<endl;
m<<"X = [ ";
for(i=0; i<nx; i++){
m<<X[i]<<" ";
}
1,2,3,4,5,6,7,8 10,11,12,13,14,15,16,17,18,19,...52
Powered by FlippingBook