NM_Lab8

#include <iostream>
#include "math.h"
#include <fstream>
#define H 0.001
using namespace std;

double a0=1.3;
double b0=1.6;
double a1y,a2y;
double c_;
double k_;
ofstream fout;

double y(double x){
    return -0.538*x*x+0.598*x+0.531+1/(-2.723*x+0.876);
}

double dy(double x){
    return -1.076*x+0.367242/((0.321704-x)*(0.321704-x))+0.598;
}

double ddy(double x){
    return 14.829458/pow((-2.723*x+0.876),3)-1.076;
}

double p(double x){
    return 1.5;
}
double q(double x){
    return -x;
}

double f(double x){
    return ddy(x)+p(x)*dy(x)+q(x)*y(x);
}

double getAij(int i, int j){
    int n=(int)((b0-a0)/H)+1;
    if(i==0){
        if(j==0)
            return c_;

        if(j==1)
            return 1;
        else
            return 0;
    }
    if(i==n-1){
        if(j==n-2)
            return 1;
        if(j==n-1)
            return k_;
        else
            return 0;
    }
    else{
        if(j==(i-1))
            return 1/(H*H)-p(a0+i*H)/(2*H);
        if(j==i)
            return q(a0+i*H)-2/(H*H);
        if(j==(i+1))
            return 1/(H*H)+p(a0+i*H)/(2*H);
        else
            return 0;
    }
}

void progonka(){
    double e=0;
    int n=(int)((b0-a0)/H)+1;
    cout << n << endl;
    double a[n];// Ay=b
    double b[n];
    double y_[n];
    a[0]=-getAij(0,1)/getAij(0,0);
    b[0]=f(a0)/getAij(0,0);
    for(int i=1; i<n-1; i++){
        a[i]=-getAij(i,i+1)/(getAij(i,i)+getAij(i,i-1)*a[i-1]);
        b[i]=(-getAij(i,i-1)*b[i-1]+f(a0+H*i))/(getAij(i,i)+getAij(i,i-1)*a[i-1]);
    }
    b[n-1]=(-getAij(n-1,n-1-1)*b[n-1-1]+f(b0))/(getAij(n-1,n-1)+getAij(n-1,n-1-1)*a[n-1-1]);
    y_[n-1]=b[n-1];//y
    for(int i=n-2; i>=0; i--){
        y_[i]=a[i]*y_[i+1]+b[i];
    }
    cout << " X    Ysr      Yi    |Ysr-Yi|" << endl;
    for( int i=0; i<n; i++){
        if(e<fabs(y_[i]-y(a0+i*H)))
            e=fabs(y_[i]-y(a0+i*H));
        fout << a0+i*H << " " << y_[i] << " " << y(a0+i*H) << " " << fabs(y_[i]-y(a0+i*H)) << endl;
        cout << a0+i*H << " " << y_[i] << " " << y(a0+i*H) << " " << fabs(y_[i]-y(a0+i*H)) << endl;
    }
    cout << "max|ei|: " << e << endl;
    fout << "max|ei|: " << e << endl;
}

int main()
{
    cout << f(a0) << " " << y(a0) << " " << y(a0+H) << endl;
    c_=(f(a0)-y(a0+H))/y(a0);
    k_=(f(b0)-y(b0-H))/y(b0);
    int n=(int)ceil((b0-a0)/H)+1;
    a1y=y(a0)-0.5*dy(a0);
    a2y=dy(b0);
    fout.open("result.txt");
    cout << "H: " << H << endl;
    progonka();
    fout.close();
    return 0;
}