usfkad6_05.cpp

Solutions are obtained for Poissson, diffusion, or wave PDEs homogeneous or nonhomogeneous equations

          //* to Latex trans1ator */
#include <string>
//#include <stdio.h>
#include <iostream>
#include <fstream>
#include <math.h>
// include <string.h>
using namespace std;

// Global Variables
int eqtype[8];
std::string mtrx[10][10];
/* the vector eqtype holds the equation parameters ( PDE type/
dimension,coordinate system, coord1, coord2, coord3, time) the coordinates are
referenced by the row number in the BC mtrx below*/
/* The boundary conditions mtrx[8][10] will hold the low end
and the high end equation parameters; these parameters are to
be entered by the user. */

char filename[20];
ofstream outFilE;
// The name of the output file.

struct   pdeans
{
	string ans1;
	string ans2;
	string ans3;
	string ans4;
	string ans5;
	string symb1;
	string symb2;
	string KAPPA;
	string nonhomf;
	string AAA;
	string Nonans1;
};

pdeans sul(pdeans, int);
pdeans sulsul(pdeans, int);
/* sul is a function that returns the solutions; sulsul is its continuation */

pdeans KandF(pdeans, std::string, int);
//  KandF prints the string for KAPPA and the nonhomogeneity f

pdeans Green(pdeans, int);
// Appends Green;s function term

pdeans KAPPAGR(pdeans);
// Coefficients for Green

pdeans Trans(pdeans, int);
//  Appends transient terms

pdeans KAPGRTR(pdeans, int);
//  Adds time factors

void main(void)
{
pdeans A;
pdeans B;
	
	
// dummyIO is used to send strings through the cin and cout channels.
std::string dummyIO;

/* dim is integer used to store the dimension 2 or 3. */ 
int dim;

/* a, i, t, and k are multipurpose integers to use in the program */
int a, i, k, t;
/* lowhigh distinguishes the low and high ends*/
int lowhigh, lohiint, p;
std::string lohisym;
/* sterm is the string value of the term */
char sterm[20];
std::string ssterm;
/* f is the string that holds the subroutine name according to the
eqn type and boundary conds; f4 holds S, O, or I */
std::string f, f4;

/* nf is the string that holds the subroutine name for the nonhomog function */ 
std::string nf;

/* to hold the solution */ 
//std::string sym1, sym2; 
//std::string solv;

/* term is the variable that holds the number of terms and the current term number */ 
int term;

/*r holds the nonhomog row number and re holds the homog row numbers */ 
int r, re;

/* These are string arrays that accumulate the answers*/
std::string Symb1, Symb2, Eigans1, Nonans1, Eigans2, Nonans2, Eigans53, Nonans3, Eigans4, Nonans4;

/* variable is the the variable names matrix */ 
std::string variable[8];

/* BC holds the questions for the Boundary conditions */
std::string BC[8];

/* the output file variable name */ 
//ofstream outFilE; 

cout << "TruKad solver: A. D. Snider and S. Kadamani.\nVersion 6/30/05\n";
cout << "See Snider's Partial Differential Equations: Sources and Solutions \n";
cout << "(Prentice Hall) for text references and notation. \n";

printf("Type the output file name (with .tex): "); 
scanf("%s",&filename); 

/* Open the file for output */ 
outFilE.open(filename); 

/* Entering the eqtype Vector values */ 

printf("Select the Partial Differential Equation\n"); 
/* printing the selections */ 
cout << "0 - Laplace or Poisson: Laplacian Psi + f(interior) = 0 . (Note sign of f!)\n";
cout << "l - Diffusion, Time Domain: d Psi/dt = Laplacian Psi + f(interior)\n";
cout << "2 - Diffusion, s-plane: s Psi - Psi(t=0) = Laplacian Psi + F(interior)\n";
cout << "3 - Wave, Time Domain: d^2 Psi/dt^2 = Laplacian Psi + f(interior)\n";
cout << "4 - Wave, s-plane: s^2Psi - sPsi(t=o) - Psi'(t=0) = Laplacian Psi + F(interior)\n";
cout << "5 - Wave, frequency domain: - omega^2 Psi = Laplacian Psi + F(interior)\n"; 
scanf("%d",&a);
eqtype[0] = a;

/* Entering nonhomogeneity flag */ 
printf("Is the PDE homogeneous (enter 0) or nonhomogeneous (enter 1)?\n"); 
scanf("%d",&a);
eqtype[7] = a;

/* Entering the dimension */ 
printf("Enter 1, 2, or 3 for 1, 2, or 3 dimensions \n"); 
scanf("%d",&a);
eqtype[1]=a;
dim=a;

/* Entering the Coordinate System */ 
printf("Select the Coordinate System\n");
printf("0 - Rectangular\nl - Cylindrical or Polar\n2 - Spherical\n");
scanf ("%d",&a);
eqtype[2]=a;

// No polar s-plane or frequency domain analyses!
/*if ((eqtype[0]==2 || eqtype[0]==4 || eqtype[0]==5) && eqtype[1]==2 && eqtype[2]==1)
{
	printf("User: the expansions for this analysis are very complicated. \n");
	printf("Please start over (CRLT c) and use a 3-dimensional analysis with\n");
	printf("homogeneous Neumann conditions on low and high z ends,\n");
	printf("and interpret the expansions with k_z equal zero only.\n");
	printf("************************************************************\n");
	printf("************************************************************\n");
    printf("\n");
    printf("\n");
}*/

/* Determining the "time" Coordinate */ 
if (eqtype[0]==1 || eqtype[0]==3)
// time domain
{
 eqtype[6]=7;
}
else
{
 /* when entry #6 is 100 in the eqtype vector it means don't care*/ 
 eqtype[6]=100;
}

/* Determining the coordinates */ 
if (eqtype[2]==0)
    {
      eqtype[3]=0;
      eqtype[4]=1;
      eqtype[5]=2;
     }
if (eqtype[2]==1)
    {
      eqtype[3]=3;
      eqtype[4]=2;
      eqtype[5]=5;
          if (eqtype[1]==2)
             {
                eqtype[4]=4;
              }     //The Cauchy-Euler equation only applies in steady state situations; this is handled later.
    }
if (eqtype[2]==2)
    {
      eqtype[5]=6;
	  eqtype[3]=8;
	  eqtype[4]=9;
    }


/* Filling the variable matrix to be used when filling the conditions matrix */
variable[0]="x"; 
variable[1]="y"; 
variable[2]="z"; 
variable[3]="\\theta"; 
variable[4]="r";
variable[5]="r";
variable[6]="r"; 
variable[7]="t";
variable[8]="\\phi";
variable[9]="\\theta";

/* The BC queries */
BC[0]= "  Select the boundary condition at the lower end for the coordinate "; 
BC[1]= "  Select the boundary condition at the upper end for the coordinate "; 
BC[2] = "\nEnter 1 for Dirichlet, Homogeneous\nEnter 2 for Dirichlet, Nonhomogeneous\nEnter 3 for Neuman, Homogeneous\n";
BC[3]="Enter 4 for Neuman, Nonhomogeneous\nEnter 5 for Robin, Homogeneous\nEnter 6 for Robin, Nonhomogeneous\n";
BC[4] = "Enter 7 for Periodic Boundry Conditions (Enter 7 for both endpoints)\n";
BC[5]="Enter 8 for Singular Boundary Condition\n";
BC[6] = "Enter 9 for Singular, Sommerfeld Outgoing Wave Condition\n";
BC[7] = "Enter 10 for Singular, Sommerfeld Incoming Wave Condition\n";
/* Requesting the BC Data */
/* Make k list the coordinates rows, 3 to 3 + #dim -1 */
//printf ( "Answer the Next questions with Upper Case Only!\n");

for (k=3;k<(3+dim);k++)
  {
   /* Filling the Lower End Boundary Conditions*/
// Jump to rho for spherical coordinates
	if(eqtype[2]==2) 
	{
		k=5;
	}

	dummyIO= BC[0]+ variable[eqtype[k]];
	cout << dummyIO << BC[2] << BC[3] << BC[4] << BC[5] << BC[6] << BC[7] << endl;
//	printf("%s ",BC[0]+ variable[eqtype[k]]);
	cin>> a;    
	if(a==1)
	{
		mtrx[eqtype[k]][0]="R";
		mtrx[eqtype[k]][1]="D";
		mtrx[eqtype[k]][2]="H";

	}
	else if(a==2)
	{
		mtrx[eqtype[k]][0]="R";
		mtrx[eqtype[k]][1]="D";
		mtrx[eqtype[k]][2]="N";
	}
	else if(a==3)
	{
		mtrx[eqtype[k]][0]="R";
		mtrx[eqtype[k]][1]="N";
		mtrx[eqtype[k]][2]="H";
	}
	else if(a==4)
	{
		mtrx[eqtype[k]][0]="R";
		mtrx[eqtype[k]][1]="N";
		mtrx[eqtype[k]][2]="N";
	}
	else if(a==5)
	{
		mtrx[eqtype[k]][0]="R";
		mtrx[eqtype[k]][1]="R";
		mtrx[eqtype[k]][2]="H";
	}
	else if(a==6)
	{
		mtrx[eqtype[k]][0]="R";
		mtrx[eqtype[k]][1]="R";
		mtrx[eqtype[k]][2]="N";
	}
	else if(a==7)
	{
		mtrx[eqtype[k]][0]="-";
		mtrx[eqtype[k]][1]="P";
		mtrx[eqtype[k]][2]="-";
	}
	else if(a==8)
	{
		mtrx[eqtype[k]][0]="-";
		mtrx[eqtype[k]][1]="S";
		mtrx[eqtype[k]][2]="-";
	}
	else if(a==9)
	{
		mtrx[eqtype[k]][0]="-";
		mtrx[eqtype[k]][1]="O";
		mtrx[eqtype[k]][2]="-";
	}
	else if(a==10)
	{
		mtrx[eqtype[k]][0]="-";
		mtrx[eqtype[k]][1]="I";
		mtrx[eqtype[k]][2]="-";
	}
	else
	{
		cout << "Invalid entry; start over" << endl;
		return;
	}

	// Filling in the upper end boundary conditions 
	dummyIO= BC[1]+ variable[eqtype[k]];
	cout << dummyIO << BC[2] << BC[3] << BC[4] << BC[5] << BC[6] << BC[7] << endl;
	cin>> a;    
	if(a==1)
	{
		mtrx[eqtype[k]][4]="R";
		mtrx[eqtype[k]][5]="D";
		mtrx[eqtype[k]][6]="H";
	}
	else if(a==2)
	{
		mtrx[eqtype[k]][4]="R";
		mtrx[eqtype[k]][5]="D";
		mtrx[eqtype[k]][6]="N";
	}
	else if(a==3)
	{
		mtrx[eqtype[k]][4]="R";
		mtrx[eqtype[k]][5]="N";
		mtrx[eqtype[k]][6]="H";
	}
	else if(a==4)
	{
		mtrx[eqtype[k]][4]="R";
		mtrx[eqtype[k]][5]="N";
		mtrx[eqtype[k]][6]="N";
	}
	else if(a==5)
	{
		mtrx[eqtype[k]][4]="R";
		mtrx[eqtype[k]][5]="R";
		mtrx[eqtype[k]][6]="H";
	}
	else if(a==6)
	{
		mtrx[eqtype[k]][4]="R";
		mtrx[eqtype[k]][5]="R";
		mtrx[eqtype[k]][6]="N";
	}
	else if(a==7)
	{
		mtrx[eqtype[k]][4]="-";
		mtrx[eqtype[k]][5]="P";
		mtrx[eqtype[k]][6]="-";
	}
	else if(a==8)
	{
		mtrx[eqtype[k]][4]="-";
		mtrx[eqtype[k]][5]="S";
		mtrx[eqtype[k]][6]="-";
	}
	else if(a==9)
	{
		mtrx[eqtype[k]][4]="-";
		mtrx[eqtype[k]][5]="O";
		mtrx[eqtype[k]][6]="-";
	}
	else if(a==10)
	{
		mtrx[eqtype[k]][4]="-";
		mtrx[eqtype[k]][5]="I";
		mtrx[eqtype[k]][6]="-";
	}
	else
	{
		cout << "Invalid entry; start over" << endl;
		return;
	}
	} // End k loop


/* Filling in BC column 3 */
  
  mtrx[0][3]="CC";
  mtrx[1][3]="CC";
  mtrx[2][3]="CC";
  mtrx[3][3]="CC";
  mtrx[4][3]="CE";
  mtrx[5][3]="BB";
  mtrx[6][3]="SB";
  mtrx[7][3]="CC";
  mtrx[8][3]="LG";
  mtrx[9][3]="SH";

// Filling in BC columns 7,8,9
  mtrx[0][7]="x";
  mtrx[0][8]="0";
  mtrx[0][9]="X";
  mtrx[1][7]="y";
  mtrx[1][8]="0";
  mtrx[1][9]="Y";
  mtrx[2][7]="z";
  mtrx[2][8]="0";
  mtrx[2][9]="Z";
  mtrx[3][7]="\\theta";
  mtrx[3][8]="0";
  mtrx[3][9]=" \\Theta";
  mtrx[4][7]="r";
  mtrx[4][8]="a";
  mtrx[4][9]="b";
  mtrx[5][7]="r";
  mtrx[5][8]="a";
  mtrx[5][9]="b";
  mtrx[6][7]="r";
  mtrx[6][8]="a";
  mtrx[6][9]="b";
  mtrx[7][7]="t";
  mtrx[7][8]="0";
  mtrx[7][9]="0";
  mtrx[8][7]="\\phi";
  mtrx[8][8]="0";
  mtrx[8][9]="\\pi";
  mtrx[9][7]="\\theta";
  mtrx[9][8]="0";
  mtrx[9][9]="2 \\pi";

  // Filling in the spherical harmonic rows
  mtrx[8][0]="-";
  mtrx[8][1]="S";
  mtrx[8][2]="H";
  mtrx[8][4]="-";
  mtrx[8][5]="S";
  mtrx[8][6]="H";
  mtrx[9][0]="-";
  mtrx[9][1]="P";
  mtrx[9][2]="H";
  mtrx[9][4]="-";
  mtrx[9][5]="P";
  mtrx[9][6]="H";

  // Filling in the Cauchy-Euler row BCs
 /* mtrx[4][0]=mtrx[5][0];
  mtrx[4][1]=mtrx[5][1];
  mtrx[4][2]=mtrx[5][2];
  mtrx[4][4]=mtrx[5][4];
  mtrx[4][5]=mtrx[5][5];
  mtrx[4][6]=mtrx[5][6];   */

/* Printing the matrix on the screen */ 

  cout << "Boundary Conditions: mtrx\n" << endl << endl;
  
  for (k=3;k<(3+dim);k++)
 { 
  printf("\n"); 
  for (i=0;i<10;i++) 
   {
	  cout << mtrx[eqtype[k]][i] << "   ";
   }
 	cout << endl << endl;
  }

/* starting the solution procedure */
/* Printing Latex preamble */
outFilE << "\\documentclass[12pt]{article}\n\\begin{document}\n";

// Time Domain Heading?
if(eqtype[0]==1 || eqtype[0]== 3)
{
	outFilE << "$\\Psi=\\Psi_{steady\\, state}+\\Psi_{transient} \\\\ \\\\\n\\Psi_{steady\\, state}=";
}
else 
{
	outFilE << "$\\Psi=";
}

// Count number of nonhomo terms and write first line.
if(eqtype[7]==1 || eqtype[0]==2 || eqtype[0]==4)
//  (If nonhomo PDE or Laplace domain; Green's function term is necessary.)
{
	term=1;
}
else
{
	term=0;
}

//   Count the number of nonhomo boundary conditions.
   for(t=0;t<8;t++)
   {
	if((mtrx[t][2]=="N"))    //(t != 4)&& omitted