Jacobi method demo

by Yi Zhang

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;
}
Jacobi method

Result of Jacobi method at different iteration

Advertisements