RichardsFastSlow/src/math/algo.cpp

#include "algo.hpp"

void Thomas(size_t n,double* a,double* b,double* c,double* d,double* x){
  //d[0]=3e-5;
  for(size_t i=1;i<n;++i){
    double w=a[i-1]/b[i-1];
    b[i]=b[i]-w*c[i-1];
    d[i]=d[i]-w*d[i-1];

  }
  x[n-1]=d[n-1]/b[n-1];
  for(int i=n-2;i>=0;--i){
    x[i]=(d[i]-c[i]*x[i+1])/b[i];
  }
}

double
norm2(double* u,size_t n){
  double r=0;
  for(size_t i=0;i<n;++i){
    double t=u[i];
    r+=(t*t);
  }
  return sqrt(r);
}

double
norm1(double* u,size_t n){
  double r=0;
  for(size_t i=0;i<n;++i){
    r+=abs(u[i]);
  }
  return r;
}

double
error2(double* u,double* v,size_t n){
  double r=0;
  for(size_t i=0;i<n;++i){
    double t=u[i]-v[i];
    r+=(t*t);
  }
  return sqrt(r);
}

void
clear(double* u,size_t n){
  for(size_t i=0;i<n;++i){
    u[i]=0;
  }
}

void
display(string name,double* u,size_t n){
  cout<<"------| "<<name<<" |------"<<endl;
  for(size_t i=0;i<n;++i){
    cout<<i<<" : "<<u[i]<<endl;
  }
  cout<<"--------------------------"<<endl;
}


double bump(double x,double delta_left,double left,double delta_right,double right){
  double left2=left-delta_left;
  double right2=right+delta_right;
  if(x>left){
    if(x<right) return 1;
    if(x>right2) return 0;
    return (right2-x)/delta_right;
  }
  else if(x<left2) return 0;
  return (x-left2)/delta_left;
}

double interpol(double x0,double x1,double t){
  return x0*(1-t)+x1*t;
}