86 lines
1.8 KiB
C++
86 lines
1.8 KiB
C++
|
|
#include <iostream>
|
|
#include <sstream>
|
|
#include <vector>
|
|
#include <cmath>
|
|
#include <lapacke.h>
|
|
using namespace std;
|
|
|
|
int n = 10;
|
|
double threshold = 1.0/sqrt(2);
|
|
double lambda(double x)
|
|
{
|
|
if(x < threshold)
|
|
{
|
|
return 1;
|
|
}
|
|
return 10;
|
|
}
|
|
|
|
double fillOffDiagonal(unsigned int i)
|
|
{
|
|
|
|
double x_l = i/double(n);
|
|
double x_u = (i+1)/double(n);
|
|
|
|
if(x_l < threshold && x_u > threshold)
|
|
{
|
|
return -n*n*(threshold-x_l+10*(x_u-threshold));
|
|
}
|
|
return -n*lambda(x_l);
|
|
}
|
|
|
|
double fillDiagonal(int i)
|
|
{
|
|
double x_l = (i-1)/double(n);
|
|
double x_u = (i+1)/double(n);
|
|
if(x_l < threshold && x_u > threshold)
|
|
{
|
|
return n*n*(threshold-x_l+10*(x_u-threshold));
|
|
}
|
|
return 2*n*lambda(x_l);
|
|
}
|
|
|
|
int main()
|
|
{
|
|
//Ignored the first Row of K cause we now that u(0)=0
|
|
|
|
int rhs = 1;
|
|
unsigned int nodes = n+1;
|
|
|
|
|
|
//Lapacke overwrites upper/lower diagonal so store twice
|
|
vector<double> upperdiagonal(n,0.0);
|
|
vector<double> lowerdiagonal(n,0.0);
|
|
vector<double> diagonal(nodes,1.0);
|
|
vector<double> F(nodes,0.0);
|
|
|
|
for(int i=0; i < n; i++)
|
|
{
|
|
upperdiagonal[i] = fillOffDiagonal(i);
|
|
lowerdiagonal[i] = fillOffDiagonal(i);
|
|
diagonal[i+1] = fillDiagonal(i+1);
|
|
}
|
|
|
|
upperdiagonal[0]=0.0;
|
|
lowerdiagonal[n-1]=0.0;
|
|
diagonal[n] = 1.0;
|
|
F[n] = 1.0;
|
|
|
|
|
|
//Tridiagonal
|
|
LAPACKE_dgtsv(LAPACK_COL_MAJOR,nodes,rhs,lowerdiagonal.data(),diagonal.data(),upperdiagonal.data(),F.data(),nodes);
|
|
|
|
|
|
cout << "Solution u_h at nodes x_0..x_" << n << ":\n";
|
|
|
|
for (unsigned int i = 0; i < nodes; ++i) {
|
|
|
|
double xi = double(i) /double(n);
|
|
cout << "u_h(" << xi << ") = " << F[i] << "\n";
|
|
// F will be overwritten with sol
|
|
}
|
|
return 0;
|
|
|
|
|
|
}
|