class px
{
   public double func(double x)
   {return 1.0;}	
}

class qx
{
   public double func(double x)
   {return Math.PI*Math.PI;}	
}

class rx
{
   public double func(double x)
   {
   return 2*Math.PI*Math.PI*Math.sin(Math.PI*x); 	   
   }	
}

// To approximate the solution to the boundary-value problem
//
//     -d(p(x)*dy/dx)/dx+q(x)y=f(x) 0<=x<=1     y(0)=0  y(1)=0;
//      
// with piecewise linear function;

// phi(x)=sum(c[i]*phi[i](x))     

public class piecewise_linear_rayleigh_ritz
{
    px p;
    qx q;
    rx r;                                                                                                                 
    int N;
    double xx[];
    double c[];
    double h[];
    
    double t;
    
    public piecewise_linear_rayleigh_ritz(int NN ,px pp,qx qq,rx rr)
    {
	N=NN;
	xx=new double[N+2];    
	c=new double[N];
	h=new double[N+1];    
	t=1.0/((double)N);
	r=rr;
	q=qq;
	p=pp;
	double[] sigma=new double[N];
	xx[0]=0;
	xx[N+1]=1;
	for(int i=0;i<N+1;i++)    
	{xx[i]=0+t*i;
	//System.out.println(xx[i]);
	}    
	
	for(int i=0;i<=N-1;i++)
	{
		h[i]=xx[i+1]-xx[i];
		//System.out.println(h[i]);
	}
}                                                                                                                                     
    
	
	public double[] sigmaf(double x)
	{
		double y[]=new double[N];
		int t=0;
		for(int i=1;i<N;i++)
		{
		if((x>=0)&&(x<xx[i]))
		{y[t]=0;}	
		else if((xx[i-1]<x)&&(x<=xx[i]))	
		{y[t]=(x-xx[i-1])/h[i-1];}	
	    else if((xx[i]<x)&&(x<=xx[i+1]))
		{y[t]=(xx[i+1]-x)/h[i];}
		else if((xx[i+1]<x)&&(x<=1))
		{y[t]=0;}
		t++;
		}
		return y;
	} 
	
	
     
	 // Approximated integral with Simpson-rule
	 public double[] Q1ap()
	 {
         double yy[]=new double[N-1];	 
		 for(int i=0;i<N-1;i++)
	     {
		 yy[i]=(h[i]/12.0)*(q.func(xx[i])+q.func(xx[i+1]));
		 //System.out.println(yy[i]);
		 }	 
		 return yy;
     }  
	    // Approximated integral with Simpson-rule
	 public double[] Q2ap()
	 {
         double yy[]=new double[N+1];	 
		 for(int i=1;i<=N;i++)
	     {
		 yy[i]=(h[i-1]/12.0)*(  3.0*q.func(xx[i]) + q.func(xx[i-1])  );
		 //System.out.println(yy[i]);
		 }	 
		 return yy;
     }
	 
   
     
	       // Approximated integral with Simpson-rule
	 public double[] Q3ap()
	 {
         double yy[]=new double[N+1];	 
		 for(int i=1;i<N+1;i++)
	     {
		 yy[i]=(h[i]/12.0)*(  3.0*q.func(xx[i]) + q.func(xx[i+1])  );
		 }	 
		 return yy;
     }
     
	 public double Q4i(double xi1,double xi,double x,double h)
	 {return (1/h)*(1/h)*p.func(x);}
	 
	 public  double Q4i_int(double xi1,double xi,double h, int n)
     {       //Newton-cotes integral 3 points
             //Simpson 1/3 rule
	    
        double hh=(xi-xi1)/n;
        double h1=hh/2.0;
        double sum=0;
        int p,i;
        double x[]=new double[3];
        double f[]=new double[3];
        for(p=0;p<n;p++)
        {  {for(i=0;i<3;i++) {x[i]=xi1+i*h1+p*hh;f[i]=Q4i(xi1,xi,x[i],h);}}
           sum+=h*(f[0]+4.0*f[1]+f[2])/6.0; 
        }
        return sum;
      }
      
      public  double newton(double xi1,double xii, double h,int n)
{       //Newton-cotes integral 2 points
        //trapezoidal rule	    
        double hh=(xii-xi1)/n;
        double h1=hh;
        double sum=0;
        int p,i;
        int nn=2*n;
        double x[]=new double[2];
        double f[]=new double[2];
        double xi[]=new double[nn];
        double fi[]=new double[nn];
        int j=0;
        for(p=0;p<n;p++)
        {  {for(i=0;i<2;i++) {x[i]=xi1+i*h1+p*h1;xi[j]=x[i];f[i]=Q4i(xi1,xii,x[i],h);fi[j]=f[i];j++;}}
           sum+=h*(f[0]+f[1])/2.0; 
        }
        return sum;
}

      
      public double[] Q4ap(int n)
      {
	      double yy[]=new double[N+2];
	      int t=0;
	      for(int i=1;i<=N;i++)
	      {yy[t]=newton(xx[i-1],xx[i],h[i-1],n);t++;}
	      return yy;
	  }
         
        
         
	     // Approximated integral with Simpson-rule
	 public double[] Q4ap()
	 {
         double yy[]=new double[N+2];	 
		 for(int i=1;i<=N+1;i++)
	     {
		 yy[i]=(h[i-1]/2.0)*(p.func(xx[i])+p.func(xx[i-1]));
		 }	 
		 return yy;
     }
	    
     public double[] Q5ap()
     {
	   double yy[]=new double[N+1];	 
		 for(int i=1;i<N+1;i++)
	     {
		 yy[i]=(h[i-1]/6.0)*( 2.0*r.func(xx[i]) + r.func(xx[i-1])  );
		 }	 
		 return yy;  
	 }
	     
	 public double[] Q6ap()
     {
	   double yy[]=new double[N+1];	 
		 for(int i=0;i<N;i++)
	     {
		 yy[i]=(h[i]/6.0)*(2.0*r.func(xx[i]) + r.func(xx[i+1]));
		 
		 }	 
		 return yy;  
	 }
	      
	   public static double[] gauss(double a[][],double b[])
       {
        int n=b.length;	  
         double x[]=new double[n];
         double carpan=0;
         double toplam=0;
          
         int i,j,k;
        
         for(k=0;k<(n-1);k++)
         {
            for(i=k+1;i<n;i++)
            {  carpan=a[i][k]/a[k][k];
               a[i][k]=0;
                for(j=k+1;j<n;j++)
                {   a[i][j]-=carpan*a[k][j]; }
                    b[i]   =b[i]   -carpan*b[k];
	        }
          }
          
          x[n-1]=b[n-1]/a[n-1][n-1];
          for(i=n-2;i>=0;i--)
          {
           toplam=0;
           for(j=i+1;j<n;j++)
           {  toplam+=a[i][j]*x[j];}
              x[i]=(b[i]-toplam)/a[i][i];
            }
           return x;
         }
  
         public double[] init_A_b()
         {
	     double[][] A=new double[N][N];    
	     double b[]=new double[N];
	     double cc[]=new double[N];
	     double[] Q1;
	     double[] Q2;
	     double[] Q3;
	     double[] Q4;
	     double[] Q5;
	     double[] Q6;
	     Q1=Q1ap();
	     Q2=Q2ap();
	     Q3=Q3ap();
	     Q4=Q4ap(1);
	     Q5=Q5ap();
	     Q6=Q6ap();
	      for(int i=0;i<N;i++)
	     {A[i][i]=Q4[i]+Q4[i+1]+Q2[i]+Q3[i];}
	     for(int i=0;i<N-1;i++)
	     {A[i][i+1]=-Q4[i+1]+Q1[i];}
	     for(int i=1;i<N-1;i++)
	     {A[i][i-1]=-Q4[i]+Q1[i-1];}    
	      for(int i=1;i<N;i++)
	     {b[i]=Q5[i]+Q6[i];}   
	     System.out.println(); 
	     c=gauss(A,b);   
	     return c;        
	     }
         
	     public double[] solution()
	     {
		     double[] cc=init_A_b();
		     double sum=0;
		     double[] depo=new double[N];
		     double[] sf=new double[N];
		     int k=0;
		     for(int j=1;j<N;j++)
		     {
			     sum=0.0;
		         sf=sigmaf(xx[j+0]);
		         for(int i=0;i<N;i++)
		         {
			     sum+=cc[i]*sf[i];    
		         }
		     depo[k]=sum;
		     k++; 
	         }  
		     return depo;
		     //System.out.println(Matrix.toString(depo));
		  }
		  
		  public double exactsol(double x)
		  {
		  return Math.sin(Math.PI*x);	  
		  }
		  
		  public String toString()
		  {
		  double[] kk=solution();
		  double depo[]=new double[N];
		  for(int i=0;i<N;i++)
		  {
		  depo[i]=exactsol(xx[i+1]);	  
		  }
		  
		  double[] xx1=new double [N-1];
		  for(int i=1;i<=N-1;i++)
		  {xx1[i-1]=xx[i];}
		  Plot pp=new Plot(xx1,depo);
		  pp.addData(xx1,kk);
		  pp.plot();
		  
		  String ss="   Exact sol             Approximation              Error  \n";
		  for(int i=0;i<N;i++)
		  {
		  ss+=exactsol(xx[i+1]) + "  " + kk[i]+ " " +Math.abs(exactsol(xx[i+1])-kk[i]) +   "\n";	  
		  }	    
		  return ss;	     
		  
		  }
		  
		  
		  
		  
	 
	                                                                                               
        public static void main(String args[])
        {
	        px pxx=new px();
	        qx qxx=new qx();
	        rx rxx=new rx();
	        int N=750;
	        
	     piecewise_linear_rayleigh_ritz rt=new piecewise_linear_rayleigh_ritz(N,pxx,qxx,rxx);   
	     
	     
	     
	     System.out.println(rt.toString());
	 
	    
         }
	                                                                                                                                                                          
                                                                                                                            
                                                                                                                                                                      
                                                                                                                      
                                                                                          
                                                                                     

         
              
                   
                             
                        
                

}


