// abc.cpp : Defines the entry point for the console application.
//


#include
#include
#include
#define NUM 10

/* The procedure will return solution to x
a - vector containing lower diagonal (store from 1 to n-1)
b - contains main diagonal (store from 0 to n-1)
c - vector containing upper diagonal (store from 0 to n-2)
d - right hand side vector (from 0 to n-1)

x - sought vector of unknowns
n - system of equations is n x n - size
*/
void TridiagSolve(const double *a, const double *b, double *c, double *d, double *x, unsigned int n)
{
int i;

/* First step: Modify the coefficients. */
c[0] /= b[0]; /* Division by zero risk. */
d[0] /= b[0]; /* Division by zero would imply a singular matrix. */
for(i = 1; i < n; i++)
{
double id = (b[i] - c[i-1] * a[i]); /* Division by zero risk. */
c[i] /= id; /* Last value calculated is redundant. */
d[i] = (d[i] - d[i-1] * a[i])/id;
}

/* Second step: Back substitution. */
x[n - 1] = d[n - 1];
for(i = n - 2; i >= 0; i--)
x[i] = d[i] - c[i] * x[i + 1];
}

void thomas(double A[][NUM], double b[], double x[])
{
double beta[NUM], alfa[NUM];
int i;
alfa[0]=A[0][1]/A[0][0];
beta[0]=b[0]/A[0][0];

for(i=1; i {
alfa[i]=(A[i][i+1])/(A[i][i]-alfa[i-1]*A[i][i-1]);
beta[i]=(b[i]+A[i][i-1]*beta[i-1])/(A[i][i]-alfa[i-1]*A[i][i-1]);
}

x[NUM-1]=beta[NUM-1];

for(i=NUM-2; i>=0; --i)
{
x[i]=beta[i]+alfa[i]*x[i+1];
}


}
double p(double x)
{
return 4;
}
double dp(double x)
{
return 0;
}

double q(double x)
{
return -8;
}



int main()
{
double x0=10, xk=20;
int alfa=100, beta=400;
double a[NUM], b[NUM], c[NUM];
double h=(xk-x0)/NUM;
double x=x0, T, Tnum[NUM], Ttridiag[NUM];
double A[NUM][NUM], A2[NUM][NUM], r[NUM], r2[NUM];
int i;


/*uzupelnienie macierzy A i wektora r dla warunkow brzegowych x0*/
r[0]=h*h*q(x0)-alfa*(h*dp(x0)*0.5-p(x0));
A[0][0]=2*p(x0);
A[0][1]=h*dp(x0)*0.5+p(x0);
/*Utworzenie wektorów a,b,c*/



/*uzupelnianie macierzy A i wektora r dla układu rownan liniowych*/
for ( i = 1; i < (NUM-1); ++i)
{
x=x+h;
A[i][i-1]=-h*dp(x)*0.5+p(x);
A[i][i]=2*p(x);
A[i][i+1]=h*dp(x)*0.5+p(x);
r[i]=h*h*q(x);

}

/*uzupelnienie macierzy A i wektora r dla warunkow brzegowych xk*/

r[NUM-1]=h*h*q(xk)-beta*(-h*dp(xk)*0.5-p(xk));
A[NUM-1][NUM-2]=-h*dp(xk)*0.5+p(xk);
A[NUM-1][NUM-1]=2*p(xk);

/*---------------------A2--------------------------------------------*/
/*uzupelnienie macierzy A i wektora r dla warunkow brzegowych x0*/
r2[0]=h*h*q(x0)-alfa*(h*dp(x0)*0.5-p(x0));
A2[0][0]=2*p(x0);
A2[0][1]=-h*dp(x0)*0.5-p(x0);
/*Utworzenie wektorów a,b,c*/



/*uzupelnianie macierzy A i wektora r dla układu rownan liniowych*/
for ( i = 1; i < (NUM-1); ++i)
{
x=x+h;
A2[i][i-1]=h*dp(x)*0.5-p(x);
A2[i][i]=2*p(x);
A2[i][i+1]=-h*dp(x)*0.5-p(x);
r2[i]=h*h*q(x);

}

/*uzupelnienie macierzy A i wektora r dla warunkow brzegowych xk*/

r2[NUM-1]=h*h*q(xk)-beta*(-h*dp(xk)*0.5-p(xk));
A2[NUM-1][NUM-2]=h*dp(xk)*0.5-p(xk);
A2[NUM-1][NUM-1]=2*p(xk);
/*------------------------------------------------------------*/
for(i=1; i {
a[i]=A2[i][i-1];
}
for(i=0; i {
b[i]=A2[i][i];
}
for(i=0; i<(NUM-1); i++)
{
c[i]=A2[i][i+1];
}


/*Rozwiazanie ukladu rownan*/
thomas(A, r, Tnum);
TridiagSolve(a, b, c, r2, Ttridiag, NUM);


/*Wyświetlanie wyników*/
printf(" X Tdokl. Tnumer. Tridiag \n");
x=x0+h;
for ( i = 0; i < NUM-1; ++i)
{
T=x*x;
printf(" %10.6lg", x);
printf(" %10.6lg", T);
printf(" %10.6lg", Tnum[i]);
printf(" %10.6lg\n", Ttridiag[i]);

x=x+h;
}

printf("Koniec\n");
system("PAUSE");
return 0;



}