[ create a new paste ] login | about

Link: http://codepad.org/ni6fysTP    [ raw code | output | fork ]

C++, pasted on Sep 26:
#include <stdio.h>
#include <iostream>
#include <math.h>
#include <cmath>
#define  N 30000
#define n 5
//#define a -1
//#define b 2

using namespace std;

double norma1;
double C,sigma0,mju,u_0,epsilon;
double R;
double *dpdR = new double [N]; 
double *HR = new double [N], *HRx = new double [N], *HRxx = new double [N], *HRxxx = new double [N];
double *PRx = new double [N];
double dR1, dR2, dR3;
double a=-10,b=20,c,d,H,THETA_;
double RO,h_0,A;
double dzetta,ro,g;
double dR11, dR22;
double *THETA_R = new double [N], *THETA_Rx = new double [N], *THETA_Rxx = new double [N];
double *THETA_Z = new double [N], *THETA_Zx = new double [N], *THETA_ZRxx = new double [N];
double *first = new double[N];
double *second = new double[N];
double *third = new double[N], *otnoshenie = new double [N];
double J, Mn, sigmaT, deltaT;
double E, J_0, r_0;
double *dHdt = new double [N];
double *Hx = new double [N], *THETA_x = new double [N];
double Ht[n][N];
double y,fx;
int j;
double norma;
double epsilon_t = 0.05;
double *raznost=new double[N*2];
double kvadrat;
double haf, Tv,pi, r, rho, Taf, Ham; 


//double f(double x)
//{
//	if (x<=0.9) 
//	{
//		return 1;
///	}
///	else return haf;
//}

double f (double x)
{
	double exp=2.7182818284;

	return 1./(1.+pow(exp,2.*1000.*x));
}

double THETA(double x)
{
	return 1.;
}



int main()
{

	c=2.;
	d=4.;
	epsilon=0.002;
	sigma0=2.;
	mju=3.;
	u_0=4.;
	C=(pow(epsilon,3.) * sigma0) /(mju*u_0);
	R=2.;
	A=3.;
	h_0=7.;
	RO=(A*h_0*epsilon)/(pow(h_0,3.)*mju*u_0);
	ro=8.;
	g=9.;
	dzetta=(ro*g*pow(h_0,2.)*epsilon) / (mju*u_0);
	J=0.5;
	sigmaT=3.;
	deltaT=4.;
	Mn=(sigmaT*deltaT*epsilon)/(mju*u_0);
	J_0=1.;
	r_0=2.;
	E=(J_0*r_0)/(ro*u_0*h_0);

	Ham=1.;
	Tv=2.;
	pi=3.14;
	r=3.;
	rho=ro;
	Taf=1.;

	haf=pow((-Ham*Tv/(6.*pi*r*rho*(Taf-Tv))),0.3333333333);
	


	H=(b-a)/n;
	THETA_=(c-d)/n;

	double delta_t=0.0002*pow(H,4.);	

	do
	{	

		for(int i=0;i<=n;i++)
		{ 
			THETA_x[i]=c+THETA_*i;
			THETA_R[i]=THETA(THETA_x[i]);
 			Hx[i]=a+H*i;
			HR[i]=f(Hx[i])+haf;
			//cout<<"H= "<<HR[i]<<"\t THETA_= "<<THETA_R[i]<<endl;
			//printf("y=%lf\t fx=%lf\n", Hx[i],HR[i]);
		}


		//***********************************************
		//|||||||||||||||||||||||||||||||||||||||||||||||
		//-------	vichislenie dpdR	-----------------
		//***********************************************
		//c||||||||||||||||||||||||||||||||||||||||||||||

		dR1=H;
		dR2=pow(H,2.);
		dR3=pow(H,3.);



		for (int i=0;i<=n-1;i++)
		{
			HRx[i] = (HR[i+1]-HR[i])/dR1;
			HRx[n] = haf;
			//cout<<"HRx= "<<HRx[i]<<"\t HRi+1= "<<HR[i+1]<<"\t HRi= "<<HR[i]<< endl;
		}
		for(int i=1;i<=n-1;i++)
		{
			HRxx[0] = 1.0;
			HRxx[n] = haf;
			
			HRxx[i] = (HR[i+1]-2*HR[i]+HR[i-1])/(dR2);
			//cout<<"HRxx= "<<HRxx[i]<<"\t HRi+1= "<<HR[i+1]<<"\t HRi= "<<HR[i]<<"\t HRi-1= "<<HR[i-1]<<endl;
		}
		for(int i=2;i<=n-1;i++)
		{
			HRxxx[0] = 1.0;
			HRxxx[1] = (1.-HR[3]+3.*HR[2])/3.;
			HRxxx[n] = haf;

			HRxxx[i] = (HR[i+1]-3.*HR[i]+3.*HR[i-1]-HR[i-2])/(dR3);
			//cout<<"HRxxx= "<<HRxxx[i]<<"\t HRi+2= "<<HR[i+2]<<"\t HRi+1= "<<HR[i+1]<<"\t HRi-1= "<<HR[i-1]<<"\t HRi-2= "<<HR[i-2]<<endl;
		}
		for (int i=0;i<=n;i++)
		{
			PRx[i] = (-3*HRx[i])/(pow(HR[i],4.)); // dPdx=-3h'/h^4 , P =1/H^3
			//cout<<"PRx= "<<HRx[i]<<"\t HRi= "<<HR[i]<<endl;
		}
		
		for (int i=0;i<=n;i++)
		{
			dpdR[i]=-C*(1. /(THETA_R[i]*THETA_R[i]) * (HRxx[i]*THETA_R[i] + HRx[i]) + HRxxx[i]) - RO*PRx[i] + dzetta*HRx[i];
			//cout<<"dpdR= "<<dpdR[i]<<endl;
		}

		/////////////////////////////////////////////////////////////////
		////////////////vichislenie evolutcionnogo yravneniia////////////
		//////////////////////v bezrazmernom vide////////////////////////
		/////////////////////////////////////////////////////////////////

		dR11=pow(THETA_,1.);
		dR22=pow(THETA_,2.);

		for (int i=1;i<=n;i++)
		{
			THETA_Rx[i] = (THETA_R[i-1]-THETA_R[i])/(dR11);
		}
		THETA_Rx[0]=THETA_Rx[1];

		for (int i=1;i<=n-1;i++)
		{
			THETA_Rxx[i] = (THETA_R[i+1]-2.*THETA_R[i]+THETA_R[i-1])/(dR22);
		}
		THETA_Rxx[0]=THETA_Rxx[1];
		THETA_Rxx[n]=THETA_Rxx[n-1];

		for (int i=0;i<=n;i++)
		{
			THETA_Zx[i] = THETA_R[i];////
		}
		for (int i=0;i<=n;i++)
		{
			THETA_ZRxx[i] = THETA_R[i];/////
		}
		cout << "START!!!" << endl;

		for (int i=0;i<=n;i++)
		{
			cout << "i=" << i << ", " << "j=" << j << endl;
			first[i] = -Mn/2.0 *( THETA_Rxx[i]*(HR[i])*(HR[i])*THETA_R[i] + 2.0*HR[i]*THETA_Rx[i]*THETA_R[i]*HRx[i] + THETA_Rx[i]*HR[i]*HR[i] );
			second[i] = -Mn/2.0 * ( 2.0*THETA_Zx[i]*HR[i]*THETA_R[i]*HRx[i]*HRx[i] + THETA_Zx[i]*HR[i]*HR[i]*HRx[i] + HR[i]*HR[i]*THETA_R[i]*HRxx[i]*THETA_Zx[i] + HR[i]*HR[i]*THETA_R[i]*HRx[i]*THETA_ZRxx[i]);
			third[i] =  -HR[i]*HR[i]*HR[i]*dpdR[i]/3.0 - HR[i]*HR[i]*HRx[i]*dpdR[i]*THETA_R[i];
			dHdt[i] = -( E*J + 1/R*( first[i] + second[i] + third[i]));

			Ht[j+1][i] = Ht[j][i] + delta_t * dHdt[i];//metod newtona
			
			raznost[i]=(Ht[j+1][i]-Ht[j][i]);
			kvadrat+=pow(raznost[i],2.);
			norma=pow(kvadrat,0.5);
			cout<<"sdg="<<norma<<endl;
		}
		j++;
	}
	while (norma<epsilon_t);

//	FILE *const_water;
	//const_water=fopen("water","r");
	unsigned char ch;
FILE *ptr;
//ptr=fopen("my_char.txt","r");


//ch=getc(prt);

if ((ptr=fopen("my_char.txt","r"))!=NULL)
{
	ch=getc(ptr);
	while (!feof(ptr))
	{
		printf("%c",ch);
		ch=getc(ptr);
		//fscanf(f, "%f", r[c]);
	}
	fclose(ptr);
}


else printf("\nFile not found!");


	delete[] HR;
	delete[] HRx;
	delete[] HRxx;
	delete[] HRxxx;
	delete[] THETA_R;
	delete[] THETA_Rx;
	delete[] THETA_Rxx;
	delete[] THETA_Zx;
	delete[] THETA_ZRxx;
	delete[] first;
	delete[] second;
	delete[] third;
	return 0;
}


Output:
START!!!
i=0, j=0
sdg=0.000981503
i=1, j=0
sdg=0.00511792
i=2, j=0
sdg=0.00524708
i=3, j=0
sdg=0.00537316
i=4, j=0
sdg=0.00549635
i=5, j=0
sdg=0.00560053
START!!!
i=0, j=1
sdg=0.00568588
i=1, j=1
sdg=0.00758676
i=2, j=1
sdg=0.00767449
i=3, j=1
sdg=0.00776123
i=4, j=1
sdg=0.00784702
i=5, j=1
sdg=0.00792034
START!!!
i=0, j=2
sdg=0.00798092
i=1, j=2
sdg=0.00943
i=2, j=2
sdg=0.00950072
i=3, j=2
sdg=0.00957093
i=4, j=2
sdg=0.00964062
i=5, j=2
sdg=0.0097004
START!!!
i=0, j=3
sdg=0.00974993
i=1, j=3
sdg=0.0109677
i=2, j=3
sdg=0.0110286
i=3, j=3
sdg=0.0110891
i=4, j=3
sdg=0.0111493
i=5, j=3
sdg=0.0112011
START!!!
i=0, j=4
sdg=0.011244
i=1, j=4
sdg=0.0123149
i=2, j=4
sdg=0.0123691
i=3, j=-1483517220

Segmentation fault


Create a new paste based on this one


Comments: