Странная работа прогонки - C++ - Ответ 5257222

#include <iostream>
#include <math.h>
#include <fstream>
#include <iomanip>
using namespace std;
 
//f(x,t)
double fin(double x, double t){
    return -(6.0+cos(4.0*x)+16.0*t*cos(4.0*x));
}
 
//mu1(t)
double muleft(double t){
    return 0.0;
}
 
//mu2(t)
double muright(double t){
    return 6.0-t*cos(4.0);
}
 
//u(x,0)
double ubegin(double x){
    return 3.0*x*x+3.0;
}
 
//Истинное решение.
double ureal(double x, double t){
    return 3.0*x*x+3.0-t*cos(4.0*x);
}
 
//s=tau/h*h - число Куранта.
int main(){
 
    const int Nh=100; const double h=1.0/double(Nh); const double s=10.0; const double tau=s*h*h; const double sigma=0.5; const int Tp=2.0/tau;
    double u1[Nh+2],u0[Nh+2],F[Nh+2],G[Nh+2],a[Nh+2],b[Nh+2],c[Nh+2],d[Nh+2],epsilon;
    ofstream res;
 
    //Нулевой слой.
    for(int i=0;i<=Nh;i++){
        double ku = double(i)*h;
 
        u0[i]=ubegin(ku);
    }
 
    cout<<endl;
 
    for(int j=1;j<=Tp;j++){
        double tu = double(j)*tau;
 
        //Прогонка начинается. Здесь a[i]u1[i-1]+b[i]u1[i]+c[i]u1[i+1]=d[i], u1[i-1]=u1[i]*F[i]+G[i].
 
        for(int k=1;k<Nh;k++){
            double ku = double(k)*h;
 
            a[k]=s*sigma;
            b[k]=-1.0-2.0*sigma*s;
            c[k]=s*sigma;
            d[k]=(u0[k-1]+u0[k+1])*s*(sigma-1.0)+u0[k]*(-1.0+2.0*s*(1.0-sigma))-tau*(sigma*fin(ku,tu)+(1.0-sigma)*fin(ku,tu-tau));  
        }
 
        //Из аппроксимации левого граничного условия следует:
        F[0]=2.0*s/(1.0+2.0*s);
        G[0]=(-2.0*s*h*muleft(tu)+u0[0]+h*h*s*fin(0.0,tu))/(2.0*s+1.0);
 
        //Искусственным образом я коэффициенты a[1] и c[n-1] обнуляю, подставляя нужные формулы для аппроксимации граничных условий, которые точно верные.
 
        d[1]-=a[1]*G[0];
        b[1]+=a[1]*F[0];
        a[1]=0.0;
 
        u1[Nh]=muright(tu);
        d[Nh-1]-=c[Nh-1]*u1[Nh];
        c[Nh-1]=0.0;
 
        for(int k=1;k<Nh;k++){
            F[k]=-c[k]/(a[k]*F[k-1]+b[k]);
            G[k]=(d[k]-a[k]*G[k-1])/(a[k]*F[k-1]+b[k]);
        }
 
        for(int l=Nh-1;l>0;l--){
            u1[l]=u1[l+1]*F[l]+G[l];
        }
 
        u1[0]=u1[1]*F[0]+G[0];
        //u1[0]=u1[1]*2.0*s/(1.0+2.0*s)+(-2.0*s*h*muleft(tu)+u0[0]+h*h*s*fin(0.0,tu))/(2.0*s+1.0);
 
        //Вывод информации о слое на времени t=2.0 и максимальной погрешности.
 
        if(j==Tp){
            epsilon=0.0;
 
            for(int s=0;s<=Nh;s++){
                double neu;
 
                neu=abs(u1[s]-ureal(double(s)*h,double(Tp)*tau));
                if(neu>epsilon) epsilon=neu;
            }
 
            res.open("C:\\Heat equation (Sigma-weight=0.5_h=100)_not mine.txt",ios::out);
            res<<"x  U0(x,t)  U1(x,t)"<<endl;
            res<<endl;
 
            for(int s=0;s<=Nh;s++){
                res<<setprecision(10)<<fixed<<h*s;
                res<<"  ";
 
                res<<setprecision(10)<<fixed<<ureal(double(s)*h,Tp*tau);
                res<<"  ";
 
                res<<setprecision(10)<<fixed<<u1[s];
                res<<endl;
            }
 
            res<<endl;
            res<<Nh<<" "<<h<<" "<<s<<" "<<tau<<" "<<Tp<<endl;
            res<<"time=";
            res<<setprecision(10)<<fixed<<Tp*tau;
            res<<", epsilon=";
            res<<setprecision(10)<<fixed<<epsilon;
    
            res.close();
        }
 
        for(int k=0;k<=Nh;k++){
            u0[k]=u1[k];
        }
    }
 
    cout<<"This is the end!"<<endl;
    return 0;
}