fork download 
  1. #define _CRT_SECURE_NO_WARNINGS
  2. # include <cstdlib>
  3. # include <iostream>
  4. # include <iomanip>
  5. #include <string>
  6.  
  7. # include <fstream>
  8. # include <ctime>
  9. # include <cmath>
  10.  
  11.  
  12. using namespace std;
  13.  
  14. int main();
  15. void dtable_data_write(ofstream &output, int m, int n, double table[]);
  16. void dtable_write(string output_filename, int m, int n, double table[],
  17. 	bool header);
  18. void f(double a, double b, double t0, double t, int n, double x[],
  19. 	double value[]);
  20. int r83_np_fa(int n, double a[]);
  21. double *r83_np_sl(int n, double a_lu[], double b[], int job);
  22. void timestamp();
  23. void u0(double a, double b, double t0, int n, double x[], double value[]);
  24. double ua(double a, double b, double t0, double t);
  25. double ub(double a, double b, double t0, double t);
  26.  
  27.  
  28. int main()
  29.  
  30.  
  31.  
  32. {
  33. 	double *a;
  34. 	double *b;
  35. 	double *fvec;
  36. 	bool header;
  37. 	int i;
  38. 	int info;
  39. 	int j;
  40. 	int job;
  41. 	double k;
  42. 	double *t;
  43. 	double t_delt;
  44. 	string t_file;
  45. 	double t_max;
  46. 	double t_min;
  47. 	int t_num;
  48. 	double *u;
  49. 	string u_file;
  50. 	double w;
  51. 	double *x;
  52. 	double x_delt;
  53. 	string x_file;
  54. 	double x_max;
  55. 	double x_min;
  56. 	int x_num;
  57. 	setlocale(LC_ALL, "Russian");
  58. 	timestamp();
  59. 	cout << "\n";
  60. 	cout << "FD1D_HEAT_IMPLICIT\n";
  61. 	cout << "  Язык разработки C++\n";
  62. 	cout << "\n";
  63. 	cout << "  Конечно-разностное решение\n";
  64. 	cout << "  Зависимое от времени одномерное уравнение теплопроводности\n";
  65. 	cout << "\n";
  66. 	cout << "    Ut - k * Uxx = F(x,t)\n";
  67. 	cout << "\n";
  68. 	cout << "  для пространственного интервала А<=X <=B с граничными условиями\n";
  69. 	cout << "\n";
  70. 	cout << "    U(A,t) = UA(t)\n";
  71. 	cout << "    U(B,t) = UB(t)\n";
  72. 	cout << "\n";
  73. 	cout << "  и временной интервал T0 <= T <= T1 с начальным условием\n";
  74. 	cout << "\n";
  75. 	cout << "    U(X,T0) = U0(X).\n";
  76. 	cout << "\n";
  77.  
  78.  
  79. 	k = 5.0E-07;
  80. 	//
  81. 	//  Set X values.
  82. 	//
  83. 	x_min = 0.0;
  84. 	x_max = 0.1;
  85. 	x_num = 2;
  86. 	x_delt = (x_max - x_min) / (double)(x_num - 1);
  87.  
  88. 	x = new double[x_num];
  89.  
  90. 	for (i = 0; i < x_num; i++)
  91. 	{
  92. 		x[i] = ((double)(x_num - i - 1) * x_min
  93. 			+ (double)(i)* x_max)
  94. 			/ (double)(x_num - 1);
  95. 	}
  96. 	// 
  97. 	//  Set T values.
  98. 	//
  99. 	t_min = 0.0;
  100. 	t_max = 5000.0;
  101. 	t_num = 20;
  102. 	t_delt = (t_max - t_min) / (double)(t_num - 1);
  103.  
  104. 	t = new double[t_num];
  105.  
  106. 	for (j = 0; j < t_num; j++)
  107. 	{
  108. 		t[j] = ((double)(t_num - j - 1) * t_min
  109. 			+ (double)(j)* t_max)
  110. 			/ (double)(t_num - 1);
  111. 	}
  112. 	//
  113. 	//  Set the initial data, for time T_MIN.
  114. 	//
  115. 	u = new double[x_num*t_num];
  116.  
  117. 	u0(x_min, x_max, t_min, x_num, x, u);
  118. 	//
  119.  
  120. 	//
  121. 	w = k * t_delt / x_delt / x_delt;
  122.  
  123. 	a = new double[3 * x_num];
  124.  
  125. 	a[0 + 0 * 3] = 0.0;
  126.  
  127. 	a[1 + 0 * 3] = 1.0;
  128. 	a[0 + 1 * 3] = 0.0;
  129.  
  130. 	for (i = 1; i < x_num - 1; i++)
  131. 	{
  132. 		a[2 + (i - 1) * 3] = -w;
  133. 		a[1 + i * 3] = 1.0 + 2.0 * w;
  134. 		a[0 + (i + 1) * 3] = -w;
  135. 	}
  136.  
  137. 	a[2 + (x_num - 2) * 3] = 0.0;
  138. 	a[1 + (x_num - 1) * 3] = 1.0;
  139.  
  140. 	a[2 + (x_num - 1) * 3] = 0.0;
  141. 	//
  142. 	//  Factor the matrix.
  143. 	//
  144. 	info = r83_np_fa(x_num, a);
  145.  
  146. 	b = new double[x_num];
  147. 	fvec = new double[x_num];
  148.  
  149. 	for (j = 1; j < t_num; j++)
  150. 	{
  151. 		//
  152. 		//  Set the right hand side B.
  153. 		//
  154. 		b[0] = ua(x_min, x_max, t_min, t[j]);
  155.  
  156. 		f(x_min, x_max, t_min, t[j - 1], x_num, x, fvec);
  157.  
  158. 		for (i = 1; i < x_num - 1; i++)
  159. 		{
  160. 			b[i] = u[i + (j - 1)*x_num] + t_delt * fvec[i];
  161. 		}
  162.  
  163. 		b[x_num - 1] = ub(x_min, x_max, t_min, t[j]);
  164.  
  165. 		delete[] fvec;
  166.  
  167. 		job = 0;
  168. 		fvec = r83_np_sl(x_num, a, b, job);
  169.  
  170. 		for (i = 0; i < x_num; i++)
  171. 		{
  172. 			u[i + j*x_num] = fvec[i];
  173. 		}
  174. 	}
  175.  
  176. 	x_file = "x.txt";
  177. 	header = false;
  178. 	dtable_write(x_file, 1, x_num, x, header);
  179.  
  180. 	cout << "\n";
  181. 	cout << "  X data written to \"" << x_file << "\".\n";
  182.  
  183. 	t_file = "t.txt";
  184. 	header = false;
  185. 	dtable_write(t_file, 1, t_num, t, header);
  186.  
  187. 	cout << "  T data written to \"" << t_file << "\".\n";
  188.  
  189. 	u_file = "u.txt";
  190. 	header = false;
  191. 	dtable_write(u_file, x_num, t_num, u, header);
  192.  
  193. 	cout << "  U data written to \"" << u_file << "\".\n";
  194.  
  195. 	cout << "\n";
  196. 	cout << "FD1D_HEAT_IMPLICIT\n";
  197. 	cout << "  Normal end of execution.\n";
  198. 	cout << "\n";
  199. 	timestamp();
  200.  
  201. 	delete[] a;
  202. 	delete[] b;
  203. 	delete[] fvec;
  204. 	delete[] t;
  205. 	delete[] u;
  206. 	delete[] x;
  207.  
  208. 	return 0;
  209. }
  210. //****************************************************************************80
  211.  
  212. void dtable_data_write(ofstream &output, int m, int n, double table[])
  213.  
  214.  
  215. {
  216. 	int i;
  217. 	int j;
  218.  
  219. 	for (j = 0; j < n; j++)
  220. 	{
  221. 		for (i = 0; i < m; i++)
  222. 		{
  223. 			output << setw(10) << table[i + j*m] << "  ";
  224. 		}
  225. 		output << "\n";
  226. 	}
  227.  
  228. 	return;
  229. }
  230. //****************************************************************************80
  231.  
  232. void dtable_write(string output_filename, int m, int n, double table[],
  233. 	bool header)
  234.  
  235.  
  236. {
  237. 	ofstream output;
  238.  
  239. 	output.open(output_filename.c_str());
  240.  
  241. 	if (!output)
  242. 	{
  243. 		cerr << "\n";
  244. 		cerr << "DTABLE_WRITE - Fatal error!\n";
  245. 		cerr << "  Could not open the output file.\n";
  246. 		return;
  247. 	}
  248.  
  249. 	if (header)
  250. 	{
  251. 		//  dtable_header_write ( output_filename, output, m, n );
  252. 	}
  253.  
  254. 	dtable_data_write(output, m, n, table);
  255.  
  256. 	output.close();
  257.  
  258. 	return;
  259. }
  260. //****************************************************************************80
  261.  
  262. void f(double a, double b, double t0, double t, int n, double x[],
  263. 	double value[])
  264.  
  265. {
  266. 	int i;
  267.  
  268. 	for (i = 0; i < n; i++)
  269. 	{
  270. 		value[i] = 0.0;
  271. 	}
  272. 	return;
  273. }
  274. //****************************************************************************80
  275.  
  276. int r83_np_fa(int n, double a[])
  277.  
  278.  
  279. {
  280. 	int i;
  281.  
  282. 	for (i = 1; i <= n - 1; i++)
  283. 	{
  284. 		if (a[1 + (i - 1) * 3] == 0.0)
  285. 		{
  286. 			cout << "\n";
  287. 			cout << "R83_NP_FA - Fatal error!\n";
  288. 			cout << "  Zero pivot on step " << i << "\n";
  289. 			return i;
  290. 		}
  291. 		//
  292. 		//  Store the multiplier in L.
  293. 		//
  294. 		a[2 + (i - 1) * 3] = a[2 + (i - 1) * 3] / a[1 + (i - 1) * 3];
  295. 		//
  296. 		//  Modify the diagonal entry in the next column.
  297. 		//
  298. 		a[1 + i * 3] = a[1 + i * 3] - a[2 + (i - 1) * 3] * a[0 + i * 3];
  299. 	}
  300.  
  301. 	if (a[1 + (n - 1) * 3] == 0.0)
  302. 	{
  303. 		cout << "\n";
  304. 		cout << "R83_NP_FA - Fatal error!\n";
  305. 		cout << "  Zero pivot on step " << n << "\n";
  306. 		return n;
  307. 	}
  308.  
  309. 	return 0;
  310. }
  311. //****************************************************************************80
  312.  
  313. double *r83_np_sl(int n, double a_lu[], double b[], int job)
  314.  
  315.  
  316. {
  317. 	int i;
  318. 	double *x;
  319.  
  320. 	x = new double[n];
  321.  
  322. 	for (i = 0; i < n; i++)
  323. 	{
  324. 		x[i] = b[i];
  325. 	}
  326.  
  327. 	if (job == 0)
  328. 	{
  329. 		//
  330. 		//  Solve L * Y = B.
  331. 		//
  332. 		for (i = 1; i < n; i++)
  333. 		{
  334. 			x[i] = x[i] - a_lu[2 + (i - 1) * 3] * x[i - 1];
  335. 		}
  336. 		//
  337. 		//  Solve U * X = Y.
  338. 		//
  339. 		for (i = n; 1 <= i; i--)
  340. 		{
  341. 			x[i - 1] = x[i - 1] / a_lu[1 + (i - 1) * 3];
  342. 			if (1 < i)
  343. 			{
  344. 				x[i - 2] = x[i - 2] - a_lu[0 + (i - 1) * 3] * x[i - 1];
  345. 			}
  346. 		}
  347. 	}
  348. 	else
  349. 	{
  350. 		//
  351. 		//  Solve U' * Y = B
  352. 		//
  353. 		for (i = 1; i <= n; i++)
  354. 		{
  355. 			x[i - 1] = x[i - 1] / a_lu[1 + (i - 1) * 3];
  356. 			if (i < n)
  357. 			{
  358. 				x[i] = x[i] - a_lu[0 + i * 3] * x[i - 1];
  359. 			}
  360. 		}
  361. 		//
  362. 		//  Solve L' * X = Y.
  363. 		//
  364. 		for (i = n - 1; 1 <= i; i--)
  365. 		{
  366. 			x[i - 1] = x[i - 1] - a_lu[2 + (i - 1) * 3] * x[i];
  367. 		}
  368. 	}
  369.  
  370. 	return x;
  371. }
  372. //****************************************************************************80
  373.  
  374. void timestamp()
  375.  
  376. /
  377. {
  378. # define TIME_SIZE 40
  379.  
  380. 	static char time_buffer[TIME_SIZE];
  381. 	const struct tm *tm;
  382. 	size_t len;
  383. 	time_t now;
  384.  
  385. 	now = time(NULL);
  386. 	tm = localtime(&now);
  387.  
  388. 	len = strftime(time_buffer, TIME_SIZE, "%d %B %Y %I:%M:%S %p", tm);
  389.  
  390. 	cout << time_buffer << "\n";
  391.  
  392. 	return;
  393. # undef TIME_SIZE
  394. }
  395. //****************************************************************************80
  396.  
  397. void u0(double a, double b, double t0, int n, double x[], double value[])
  398.  
  399.  
  400. {
  401. 	int i;
  402.  
  403. 	for (i = 0; i < n; i++)
  404. 	{
  405. 		value[i] = 100.0;
  406. 	}
  407. 	return;
  408. }
  409. //****************************************************************************80
  410.  
  411. double ua(double a, double b, double t0, double t)
  412.  
  413.  
  414. {
  415. 	double value;
  416.  
  417. 	value = 20.0;
  418.  
  419. 	return value;
  420. }
  421.  
  422.  
  423. double ub(double a, double b, double t0, double t)
  424.  
  425. {
  426. 	double value;
  427.  
  428. 	value = 20.0;
  429. 	system("PAUSE");
  430. 	return value;
  431.  
  432.  
  433.  
  434.  
  435. }
Compilation error #stdin compilation error #stdout 0s 0KB
comments (?)
stdin 
Standard input is empty
compilation info 
prog.cpp:376:1: error: expected initializer before ‘/’ token
 /
 ^
stdout 
Standard output is empty
https://ideone.com/2CUNoK
language:
C++ (gcc 8.3)
created:
8 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!