// 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]; }
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*/
/*---------------------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]; }