### Jacobi method demo

Code below is to solve the static version of the equation in previous FDM code demo, with B.C. changed as homogeneous. The Jacobi method is used to replace the LAPACK solver which gives inexact method. The exact solution is
$u(x)=-x^2/2+x/2$
Using initial guess as
$u_0(x)=\sin(2\pi x)+0.5\sin(10\pi x)+0.3\sin(20\pi x)+0.1\sin(40\pi x)$
one can see that most of the computing cost is spent on smoothing lowest mode.


#include <stdio.h>
#include <math.h>
#include <iostream>
#include <fstream>
//#include <vector>
//#include <sunperf.h>
using namespace std;
//typedef std::vector<float> fvector_s;

//Maximum of an array.
float max(float a[], int n) {
float v = a[0];
for (int i=1; i<n; i++) {
if (a[i] > v) {
v = a[i];
}
}
return v;
}

int main ()
{
int n=500;
float tmp, b[n], a_d[n], a_u[n-1], a_l[n-1],    \
tol=1.e-7, x0[n], x1[n], x2[n], error=1.;
float h=1./(n+1), pi=3.1416;
int count=0;

//initializing a
for(int i=0; i<n; ++i){
a_d[i]=2.0;
}
for(int i=0; i<n-1; i++){
a_l[i]=-1.;
a_u[i]=-1.;
}

//initializaing b[]
for(int i=0; i<n; ++i){
b[i]=h*h;
}

//initial guess
for(int i=0; i<n; ++i){
x0[i]=sin(2*pi*h*(i+1))+0.5*sin(10*pi*h*(i+1))+\
0.3*sin(20*pi*h*(i+1))+0.1*sin(40*pi*h*(i+1));
x1[i]=x0[i];
}

ofstream outfile;
outfile.open("jac.out");
do{
count+=1;
for(int i=0; i<n; ++i){
if(i==0)
x2[i]=(b[i]-a_u[i]*x1[i+1])/a_d[i];
else if(i==n-1)
x2[i]=(b[i]-a_l[i-1]*x1[i-1])/a_d[i];
else
x2[i]=(b[i]-a_u[i]*x1[i+1]-a_l[i-1]*x1[i-1])/a_d[i];
}
for(int i=0; i<n; i++){
x1[i]=x2[i]-x1[i];
}
error=max(x1, n)/max(x0, n);
for(int i=0; i<n; i++){
x1[i]=fabs(x2[i]-(-(i*h)*(i*h)/2+i*h/2));
}
outfile << count << "\t" << error << endl;
for(int i=0; i<n; ++i){
x1[i]=x2[i];
}
}while(error>tol);
cout << "Number of iterations:" << count << endl;
for(int i=0; i<n; i++){
outfile << count << "\t" << x2[i] << endl;
}

outfile.close();
return 0;
}


Result of Jacobi method at different iteration