 |
Постоянные коэффициенты, условия 3-го рода.
|
|
|
|
Лабораторная работа №10. Метод конечных разностей для нестационарных уравнений с постоянными коэффициентами.
Постоянные коэффициенты, условия 1-го рода.
#include<stdio.h>
#include<conio.h>
#include<math.h>
#include<iostream>
#include<locale>
#include <fstream>
using namespace std;
int N, M, i, j, l, Tt;
double h, t, y, s, mAx, tim, a[20000], b[20000], X[20000], Y[20000];
double fi (int i, int j)
{ double x, y;
x=i*h; y=t*(j+0.5);
return(3*x*x*x*x*y*y-12*x*x*y*y*y);
}
double u (double x, double y)
{ return(x*x*x*x*y*y*y);
}
double A(double s)
{ return(t*s/(h*h));
}
double B(double s)
{ return(t*s/(h*h));
}
double C(double s)
{ return(1+2*t*s/(h*h));
}
double F(int i, int j)
{ double S;
S=X[i]+t*(1-s)*(X[i+1]-2*X[i]+X[i-1])/(h*h)+t*fi(i,j);
return(S);
}
void inp(double &l, double &Tt)
{ cout << "l = "; cin >> l;
cout << "T = "; cin >> Tt;
}
void out(double s, double mAx, double tim)
{ cout << "\nN = " << N << ", M = " << M << ", max = " << mAx << ", time = " << tim << ".";
}
void yav_metod(double &mAx)
{
h=l/N; t=h*h/2; M=2*Tt*N*N/(l*l);
for (i=0; i<=N; i++) { X[i]=u(i*h,0); Y[i]=0; }
for (j=1; j<=M; j++)
{ for (i=0; i<=N; i++) X[i]=Y[i];
for (i=1; i<N; i++)
{ Y[i]=X[i]+t*(X[i+1]-2*X[i]+X[i-1])/(h*h)+t*fi(i,j-1);
if (mAx<abs(u(i*h,j*t)-Y[i])) mAx=abs(u(i*h,j*t)-Y[i]); }
Y[0]= u(0,j*t); Y[N]= u(l,j*t);
}
}
void neyav_metod(double &mAx)
{
M=N; h=l/N; t=Tt/M;
for (i=0; i<=N; i++) { Y[i]=0; X[i]=u(i*h,0);}
for (j=1; j<=M; j++)
{ a[1]=0; b[1]=u(0, j*t);
Y[0]= u(0,j*t);
for (i=2; i<=N; i++)
{ b[i]=(A(s)*b[i-1]+F(i-1,j-1))/(C(s)-a[i-1]*A(s));
a[i]=B(s)/(C(s)-a[i-1]*A(s));
}
Y[N]= u(l,j*t);
for (i=N; i>=1; i--) Y[i-1]=a[i]*Y[i]+b[i];
for (i=0; i<=N; i++)
{ if (mAx<abs(Y[i]-u(h*i, t*j))) mAx=abs(Y[i]-u(h*i, t*j));
X[i]=Y[i];
}
}
}
void main ()
{
setlocale(LC_ALL, "rus");
ofstream fout("cpp.txt");
inp(l,Tt);
for (s=0; s<=1; s+=0.1)
{ cout << "\nСигма = " << s;
fout << "\nСигма = " << s;
for (N=10; N<=1000; N=N*10)
{ tim=clock(); mAx=0;
if (s==0) yav_metod(mAx); else neyav_metod(mAx);
tim=clock()-tim; out(s,mAx,tim);
fout << "\nN = " << N << ", M = " << M << "max = " << mAx << ", time = " << tim << ".";
}
}
fout.close();
|
|
|
getch();
}
Результаты выполнения программы.
l = 1, T =1.
Сигма = 0
N = 10, M = 199, max = 0.00131501, time = 1.
N = 100, M = 20000, max = 1.34993e-005, time = 805.
N = 1000, M = 2000000, max = 1.35009e-007, time = 823232.
Сигма = 0.1
N = 10, M = 10, max = 40.8816, time = 0.
N = 100, M = 100, max = 8.78392e+083, time = 12.
N = 1000, M = 1000, max = 3.40913e+305, time = 1344.
Сигма = 0.2
N = 10, M = 10, max = 0.0521667, time = 0.
N = 100, M = 100, max = 3.23156e+049, time = 12.
N = 1000, M = 1000, max = 1.#INF, time = 1296.
Сигма = 0.3
N = 10, M = 10, max = 0.0193234, time = 1.
N = 100, M = 100, max = 2.46777e+026, time = 12.
N = 1000, M = 1000, max = 1.#INF, time = 1217.
Сигма = 0.4
N = 10, M = 10, max = 0.00708031, time = 0.
N = 100, M = 100, max = 2.22503e+007, time = 12.
N = 1000, M = 1000, max = 7.08037e+161, time = 1187.
Сигма = 0.5
N = 10, M = 10, max = 0.00491127, time = 0.
N = 100, M = 100, max = 4.93183e-005, time = 12.
N = 1000, M = 1000, max = 4.93244e-007, time = 1193.
Сигма = 0.6
N = 10, M = 10, max = 0.0162144, time = 0.
N = 100, M = 100, max = 0.00122582, time = 12.
N = 1000, M = 1000, max = 0.000118422, time = 1194.
Сигма = 0.7
N = 10, M = 10, max = 0.0270821, time = 0.
N = 100, M = 100, max = 0.00239837, time = 13.
N = 1000, M = 1000, max = 0.00023631, time = 1194.
Сигма = 0.8
N = 10, M = 10, max = 0.0375311, time = 0.
N = 100, M = 100, max = 0.00356634, time = 12.
N = 1000, M = 1000, max = 0.000354153, time = 1194.
Сигма = 0.9
N = 10, M = 10, max = 0.0475779, time = 0.
N = 100, M = 100, max = 0.00472975, time = 12.
N = 1000, M = 1000, max = 0.000471949, time = 1193.
Сигма = 1
N = 10, M = 10, max = 0.057239, time = 0.
N = 100, M = 100, max = 0.00588861, time = 13.
N = 1000, M = 1000, max = 0.0005897, time = 1195.
Переменные коэффициенты, условия 1-го рода.
#include<stdio.h>
#include<conio.h>
#include<math.h>
#include<iostream>
#include<locale>
#include <fstream>
using namespace std;
int N, M, i, j, l, Tt;
double h, t, y, s, mAx, tim, a[20000], b[20000], X[20000], Y[20000];
double fi (int i, int j)
{ double x, y;
x=i*h; y=t*(j+0.5);
return(3*x*x*x*x*y*y-16*x*x*x*y*y*y*y*y);
}
double k(int i, int j)
{ return(h*(i+0.5)*t*j*t*j);
}
double u (double x, double y)
{ return(x*x*x*x*y*y*y);
}
double A(double s, int i, int j)
{ return(t*k(i-1,j)*s/(h*h));
}
double B(double s, int i, int j)
{ return(t*k(i,j)*s/(h*h));
}
double C(double s, int i, int j)
{ return(1+t*(k(i-1,j)+k(i,j))*s/(h*h));
}
double F(int i, int j)
{ double S;
S=X[i]+t*(1-s)*(k(i+1,j)*X[i+1]-(k(i+1,j)+k(i,j))*X[i]+k(i,j)*X[i-1])/(h*h)+t*fi(i,j);
return(S);
}
void inp(double &l, double &Tt)
{ cout << "l = "; cin >> l;
cout << "T = "; cin >> Tt;
}
void out(double s, double mAx, double tim)
{ cout << "\nN = " << N << ", M = " << M << ", max = " << mAx << ", time = " << tim << ".";
|
|
|
}
void yav_metod(double &mAx)
{
h=l/N; t=h*h/2; M=2*Tt*N*N/(l*l);
for (i=0; i<=N; i++) { X[i]=u(i*h,0); Y[i]=0; }
for (j=1; j<=M; j++)
{ for (i=0; i<=N; i++) X[i]=Y[i];
for (i=1; i<N; i++)
{ Y[i]=X[i]+t*(k(i,j-1)*X[i+1]-(k(i,j-1)+k(i-1,j-1))*X[i]+k(i-1,j-1)*X[i-1])/(h*h)+t*fi(i,j-1);
if (mAx<abs(u(i*h,j*t)-Y[i])) mAx=abs(u(i*h,j*t)-Y[i]); }
Y[0]= u(0,j*t); Y[N]= u(l,j*t);
}
}
void neyav_metod(double &mAx)
{
M=N; h=l/N; t=Tt/M;
for (i=0; i<=N; i++) { Y[i]=0; X[i]=u(i*h,0);}
for (j=1; j<=M; j++)
{ a[1]=0; b[1]=u(0, j*t);
Y[0]= u(0,j*t);
for (i=2; i<=N; i++)
{ b[i]=(A(s,i,j-1)*b[i-1]+F(i-1,j-1))/(C(s,i,j-1)-a[i-1]*A(s,i,j-1));
a[i]=B(s,i,j-1)/(C(s,i,j-1)-a[i-1]*A(s,i,j-1));
}
Y[N]= u(l,j*t);
for (i=N; i>=1; i--) Y[i-1]=a[i]*Y[i]+b[i];
for (i=0; i<=N; i++)
{ if (mAx<abs(Y[i]-u(h*i, t*j))) mAx=abs(Y[i]-u(h*i, t*j));
X[i]=Y[i];
}
}
}
void main ()
{
setlocale(LC_ALL, "rus");
ofstream fout("cpp.txt");
inp(l,Tt);
for (s=0; s<=1; s+=0.1)
{ cout << "\nСигма = " << s;
fout << "\nСигма = " << s;
for (N=10; N<=1000; N=N*10)
{ tim=clock(); mAx=0;
if (s==0) yav_metod(mAx); else neyav_metod(mAx);
tim=clock()-tim; out(s,mAx,tim);
fout << "\nN = " << N << ", M = " << M << "max = " << mAx << ", time = " << tim << ".";
}
}
fout.close();
getch();
}
Результаты выполнения программы.
l = 1, T =1.
Сигма = 0
N = 10, M = 199, max = 0.00208381, time = 0.
N = 100, M = 20000, max = 2.1633e-005, time = 453.
N = 1000, M = 2000000, max = 2.16372e-007, time = 457125.
Сигма = 0.1
N = 10, M = 10, max = 0.0887729, time = 0.
N = 100, M = 100, max = 3.68694e+066, time = 16.
N = 1000, M = 1000, max = 4.82742e+305, time = 2015.
Сигма = 0.2
N = 10, M = 10, max = 0.0367842, time = 0.
N = 100, M = 100, max = 1.98845e+037, time = 0.
N = 1000, M = 1000, max = 2.54744e+305, time = 1563.
Сигма = 0.3
N = 10, M = 10, max = 0.0228462, time = 0.
N = 100, M = 100, max = 1.3498e+017, time = 15.
N = 1000, M = 1000, max = 8.44708e+304, time = 797.
Сигма = 0.4
N = 10, M = 10, max = 0.0106661, time = 0.
N = 100, M = 100, max = 0.198424, time = 0.
N = 1000, M = 1000, max = 1.51908e+139, time = 703.
Сигма = 0.5
N = 10, M = 10, max = 0.0062122, time = 0.
N = 100, M = 100, max = 0.000516836, time = 16.
N = 1000, M = 1000, max = 5.30847e-005, time = 687.
Сигма = 0.6
N = 10, M = 10, max = 0.0162701, time = 0.
N = 100, M = 100, max = 0.00144204, time = 0.
N = 1000, M = 1000, max = 0.000141908, time = 703.
Сигма = 0.7
N = 10, M = 10, max = 0.027363, time = 0.
N = 100, M = 100, max = 0.00265983, time = 16.
N = 1000, M = 1000, max = 0.000264725, time = 687.
Сигма = 0.8
N = 10, M = 10, max = 0.0382426, time = 0.
N = 100, M = 100, max = 0.00390809, time = 0.
N = 1000, M = 1000, max = 0.000391118, time = 703.
Сигма = 0.9
N = 10, M = 10, max = 0.0487114, time = 0.
N = 100, M = 100, max = 0.00516122, time = 16.
N = 1000, M = 1000, max = 0.000518473, time = 703.
Сигма = 1
N = 10, M = 10, max = 0.0587916, time = 0.
N = 100, M = 100, max = 0.00641366, time = 0.
N = 1000, M = 1000, max = 0.000646185, time = 703.
.
Постоянные коэффициенты, условия 3-го рода.
#include<stdio.h>
#include<conio.h>
#include<math.h>
#include<iostream>
#include<locale>
#include <fstream>
using namespace std;
int N, M, i, j, l, Tt;
double h, t, y, s, mAx, tim, b1, b2, x2, v2, a[20000], b[20000], X[20000], Y[20000];
|
|
|
double fi (double x, double y)
{ return(3*x*x*x*x*y*y-12*x*x*y*y*y);
}
double u (double x, double y)
{ return(x*x*x*x*y*y*y);
}
double A(double s)
{ return(t*s/(h*h));
}
double B(double s)
{ return(t*s/(h*h));
}
double C(double s)
{ return(1+2*t*s/(h*h));
}
double F(int i, int j)
{ double S;
S=X[i]+t*(1-s)*(X[i+1]-2*X[i]+X[i-1])/(h*h)+t*fi(i*h,t*(j-0.5));
return(S);
}
double u1(int j)
{ return(0.5*h*fi(0,t*(j-0.5)));
}
double u2(int j)
{ double y;
y=t*(j-0.5);
return(4*y*y*y*l*l*l+0.5*h*fi(l,y));
}
void inp(double &l, double &Tt, double &b1, double &b2)
{ cout << "l = "; cin >> l;
cout << "T = "; cin >> Tt;
b1=1; b2=0;
}
void out(double s, double mAx, double tim)
{ cout << "\nN = " << N << ", M = " << M << ", max = " << mAx << ", time = " << tim << ".";
}
void yav_metod(double &mAx)
{
h=l/N; t=h*h/2; M=2*Tt*N*N/(l*l);
for (i=0; i<=N; i++) { X[i]=u(i*h,0); Y[i]=0; }
for (j=1; j<=M; j++)
{ for (i=0; i<=N; i++) X[i]=Y[i];
for (i=1; i<N; i++)
{ Y[i]=X[i]+t*(X[i+1]-2*X[i]+X[i-1])/(h*h)+t*fi(i*h,t*(j-0.5));
if (mAx<abs(u(i*h,j*t)-Y[i])) mAx=abs(u(i*h,j*t)-Y[i]);
}
Y[0]=t*((X[1]-X[0])/h-b1*X[0]+u1(j))/(0.5*h)+X[0];
Y[N]=t*((X[N-1]-X[N])/h-b2*X[N]+u2(j))/(0.5*h)+X[N];
}
}
void neyav_metod(double &mAx)
{ M=N; h=l/N; t=Tt/M;
for (i=0; i<=N; i++) { Y[i]=0; X[i]=u(i*h,0);}
for (j=1; j<=M; j++)
{ a[1]=s/(s+s*b1*h+h*h/(2*t));
b[1]=((1-s)*X[1]/h+(h/(2*t)-(1-s)/h-(1-s)*b1)*X[0]+u1(j))/(s/h+s*b1+h/(2*t));
for (i=2; i<=N; i++)
{ b[i]=(A(s)*b[i-1]+F(i-1,j))/(C(s)-a[i-1]*A(s));
a[i]=B(s)/(C(s)-a[i-1]*A(s));
}
x2=s/(s+s*b2*h+h*h/(2*t));;
v2=((1-s)*X[N-1]/h+(h/(2*t)-(1-s)/h-(1-s)*b2)*X[N]+u2(j))/(s/h+s*b2+h/(2*t));
Y[N]=(v2+x2*b[N])/(1-a[N]*x2);
for (i=N; i>=1; i--) Y[i-1]=a[i]*Y[i]+b[i];
for (i=0; i<=N; i++)
{ if (mAx<abs(Y[i]-u(h*i, t*j))) mAx=abs(Y[i]-u(h*i, t*j));
X[i]=Y[i];
}
}
}
void main ()
{
setlocale(LC_ALL, "rus");
ofstream fout("cpp.txt");
inp(l,Tt,b1,b2);
for (s=0; s<=1; s+=0.1)
{ cout << "\nСигма = " << s;
fout << "\nСигма = " << s;
for (N=10; N<=1000; N=N*10)
{ tim=clock(); mAx=0;
if (s==0) yav_metod(mAx); else neyav_metod(mAx);
tim=clock()-tim; out(s,mAx,tim);
fout << "\nN = " << N << ", M = " << M << "max = " << mAx << ", time = " << tim << ".";
}
}
fout.close();
getch();
}
Результаты выполнения программы.
l = 1, T =1.
Сигма = 0
N = 10, M = 199, max = 0.00928976, time = 0.
N = 100, M = 20000, max = 0.000104699, time = 218.
N = 1000, M = 2000000, max = 1.05609e-006, time = 208954.
Сигма = 0.1
N = 10, M = 10, max = 3353.23, time = 0.
N = 100, M = 100, max = 1.40445e+085, time = 12.
N = 1000, M = 1000, max = 2.82242e+305, time = 1254.
Сигма = 0.2
N = 10, M = 10, max = 5.14895, time = 0.
N = 100, M = 100, max = 6.62593e+050, time = 10.
N = 1000, M = 1000, max = 3.58088e+305, time = 1062.
Сигма = 0.3
N = 10, M = 10, max = 0.0540937, time = 0.
N = 100, M = 100, max = 4.01059e+027, time = 10.
|
|
|
N = 1000, M = 1000, max = 1.#INF, time = 826.
Сигма = 0.4
N = 10, M = 10, max = 0.0071138, time = 0.
N = 100, M = 100, max = 1.8698e+008, time = 8.
N = 1000, M = 1000, max = 5.7229e+162, time = 709.
Сигма = 0.5
N = 10, M = 10, max = 0.0224186, time = 0.
N = 100, M = 100, max = 0.000226458, time = 7.
N = 1000, M = 1000, max = 2.2648e-006, time = 677.
Сигма = 0.6
N = 10, M = 10, max = 0.0438244, time = 0.
N = 100, M = 100, max = 0.00246408, time = 7.
N = 1000, M = 1000, max = 0.000226469, time = 869.
Сигма = 0.7
N = 10, M = 10, max = 0.0644295, time = 0.
N = 100, M = 100, max = 0.00469324, time = 12.
N = 1000, M = 1000, max = 0.000450589, time = 1104.
Сигма = 0.8
N = 10, M = 10, max = 0.08426, time = 0.
N = 100, M = 100, max = 0.00691396, time = 10.
N = 1000, M = 1000, max = 0.000674623, time = 997.
Сигма = 0.9
N = 10, M = 10, max = 0.103343, time = 0.
N = 100, M = 100, max = 0.00912629, time = 8.
N = 1000, M = 1000, max = 0.000898573, time = 819.
Сигма = 1
N = 10, M = 10, max = 0.121704, time = 0.
N = 100, M = 100, max = 0.0113302, time = 7.
N = 1000, M = 1000, max = 0.00112244, time = 745.
|
|
|
