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
|
|