Finite difference method for partial differential equations

Hello, I'm from Russia, need your help in finding bugs. Calculated by the explicit scheme. Produces some very large numbers.

task:
[math] U_t = 3 (1,1-0,5 x) U_ {xx} + e ^ t-1 [/ math]
[math] U (0, t) = 0 [/ math]
[math] U (1, t) = 0 [/ math]
[math] U (x, 0) = 0.01 (1-x) x [/ math]

Need to find a solution with accuracy [math] 0.0001 [/ math] on the interval [math] T = 1 / a ^ *, where a ^ * = \ max a (x, t) [/ math]
Plot graphs of functions [math] u (x ^ *, t), u (x, jt ^ *) [/ math] where [math] x ^ * = 0.6, t ^ * = T/10, j = 1,2,4 [/ math]

explicit difference scheme is as follows:
([math] $ u_t ^ {j +1}-u_i ^ j) / \ tau = 3 (1,1-0,5 x_i) (u_ {i +1} ^ {j}-2u_i ^ j + u_ {i -1} ^ j) / h ^ 2 + e ^ {t_j} +1 $ $ [/ math]

code of the program:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
int main ( void )
{
	setlocale(LC_ALL, "rus");

	int I = 10, J = 30, i, j;
	double  T = 1.0/ pow(3.3, 0.5), h_x = 1.0/ I, h_t = T/ J, epsilon = h_t + pow(h_x, 2), c;
	double **u = new double *[I + 1];
	for (i = 0; i <= I; i++) u[i] = new double [J + 1]; 

	cout<< "Схема может быть неустойчива  при значениях Х :\n";
	for (i = 0; i <= I; i++)
	{
		c = 3 * (1.1 - 0.5 * h_x * i) * h_t * pow(h_x, -2);
		if (c < 0.5) cout << i * h_x << "   ";
	}
	cout <<"\n";
	
	//нулевой слой (j = 0)
 	for (i = 0; i <= I; i++)
	{
		u [i][0] = 0.01 * (1 - i * h_x) * i * h_x;
		//u [i][0] = 1 - i * h_x;  //НУ, несоответствие ГУ и НУ!
	}
		//последующие слои
	for (j = 0; j <= J; j++)
	{
		for (i = 1; i < I; i++) //расчёт j + 1 - го слоя по j-му
		{
			u [0][j + 1] = 0; //ГУ u [0][j + 1] = 1;
			u [I][j + 1] = 0; //ГУ
			u [i][j + 1] = u [i][j] + h_t * (3 * (1.1 - 0.5 * h_x * i) * (u [i + 1][j] -2 * u [i][j] + u [i - 1][j])/ pow(h_x, 2) + exp(h_t * j) - 1);
		}
	}
	int Jv = J/10;

	ofstream out;
    out.open ("D:\\proga7.txt");

	out << "U = U(0.6, t):\n";
	cout << "U = U(0.6, t):\n";
	for (i = 0; i <= J; i++)
	{
		out << h_t * i <<"\t"<< u [6][i] <<"\n";
	    cout << h_t * i <<"\t"<< u [6][i] <<"\n";
	}
	out << "\n U = U(x, 0.33):\n";
	cout << "\n U = U(x, 0.33):\n";
	for (i = 0; i <= I; i++)
	{
		out << h_x * i <<"\t"<< u [i][Jv] <<"\n";
		cout << h_x * i <<"\t"<< u [i][Jv] <<"\n";
	}
	out << "\n U = U(x, 0.66):\n";
	cout << "\n U = U(x, 0.66):\n";
	for (i = 0; i <= I; i++)
	{
		out << h_x * i <<"\t"<< u [i][Jv * 2] <<"\n";
	    cout << h_x * i <<"\t"<< u [i][Jv * 2] <<"\n";
	}
	out << "\n U = U(x, 1.32):\n";
	cout << "\n U = U(x, 1.32):\n";
	for (i = 0; i <= I; i++)
	{
		out << h_x * i <<"\t"<< u [i][Jv * 4] <<"\n";
	    cout << h_x * i <<"\t"<< u [i][Jv * 4] <<"\n";
	}
	out.close();
	getch();
	return 0;
}


displays the following:
[url=http://imageshack.us/photo/my-images/585/s8mi.png/][img]http://img585.imageshack.us/img585/4125/s8mi.png[/img][/url][/quote]
Use a smaller time step.
¿what's the criteria that you are using in line 14?

Also line 25 should be for (j = 0; j <J; j++) or you will be accessing out of bounds (you try to access j+1)


Some complains:
_ When posting code try to make it self-contained. You are missing headers.
_ Limit the scope of your variables. `i' and `j' are only used for traversing the matrix, they should only exist inside the loops.
_ You are leaking memory. There is no reason to use dynamic allocation here, but if you must (when you increase the size of the matrix), use a container like std::vector or std::valarray.
_ I, J, T, c are not good variable names.
_ Use whitespace to group, by instance (u[i+1][j] - 2*u[i][j] + u[i-1][j])
_ No need to use `pow()' to find an square h_x*h_x
_ It's bothersome that you put your output as an image. It is text, so put it as text. Use [output][/output] tags
спасибо!
Topic archived. No new replies allowed.