Численное решение краевых задач математической физики методом сеток

Рисунок 2 – График решения при фиксированной координате

Рисунок 3 – График решения при фиксированной координате

Рисунок 4 – График решения при фиксированном времени

Рисунок 5 – График решения при фиксированном времени

7 Экспериментальное исследование скорости сходимости

Для экспериментального исследования скорости сходимости явной схемы  изучим скорость изменения погрешности  при уменьшении шагов сетки.

Зададим погрешность сеточного  решения относительно аналитического следующим образом:

где и - некоторые константы, а .

Спланируем эксперимент  следующим образом:

Зафиксируем определенное количество шагов сетки (по координате или по времени), а по другой оси его будем увеличивать. На каждой итерации процесса будем вычислять получившуюся погрешность. Погрешность будем рассчитывать как разницу с аналитическим решением в заданной точке, посчитанной с точностью

Результаты представим на графиках, чтобы визуально оценить характер сходимости, а также в таблицах.

Зафиксируем 400 разбиений  по пространственному измерению, получим

Рисунок  6– Зависимость погрешности решения от мелкости сетки по

_{Приведем  таблицу значений погрешностей для  этого случая}

K	I
10	400	3.0	0.058
20	400	1.5	0.041	1.41
40	400	0.75	0.029	1.41
80	400	0.375	0.0205	1.41
160	400	0.188	0.0146	1.40
320	400	0.094	0.0104	1.36
640	400	0.047	0.0075	1.389
1280	400	0.023	0.00557	1.35
2560	400	0.012	0.0042	1.308
5120	400	0.006	0.00335	1.272

Таблица 1  – исследование скорости сходимости.

Для анализа зависимости  от величины шага по пространственной координате, будем рассчитывать количество шагов по временному измерению по эвристической формуле чтобы гарантировать, что погрешность, получаемая за счет величины шага по времени, заведомо меньше, чем результирующая. Приведем график и таблицу результатов.

Рисунок  7– Зависимость погрешности решения от мелкости сетки по

 

K	I
600	10	0.3	0.249
2100	20	0.15	0.125	1.992
8100	40	0.075	0.062	2.01
32100	80	0.038	0.028	2.214

Таблица 2 – исследование скорости сходимости.

 

По итогам эксперимента можно  сделать вывод о справедливости нашей аналитической оценки скорости сходимости от величины шага по измерению координаты, но с зависимостью от величины шага по времени не все очевидно. По данным в таблице видно, что скорость сходимости не достигает полученных аналитически оценок, и, мало того, убывает при измельчении сетки.

 

Заключение

В процессе выполнения работы методом конечных разностей была построена разностная схема Кранка-Николсона для исходной задачи математической физики. Для схемы было проведено исследование аппроксимации, была установлена условная устойчивость и выполнение необходимого признака Неймана.

Был создан программный модуль, моделирующий процесс распространения тепла в тонком стержне. Для проверки правильности подсчета порядка аппроксимации был проведён вычислительный эксперимент с целью установить зависимость погрешности численного решения задачи от шага сетки. При моделировании была подтверждена сходимость решения к истинному на рассмотренных сетках, однако остался незакрытым вопрос о скорости этой сходимости при измельчении шага по времени, т.к. при эксперименте аналитическая оценка не подтвердилась.

В результате были выделены следующие достоинства метода конечных разностей:

• широкая распространенность (большое количество доступных источников информации);
• сходимость к аналитическому решению (в данной работе она наблюдалась);
• подходит для решения класса задач, исследованного в данной работе;
• простота построения дискретных моделей.

Отметим, что конечных разностей не лишён следующих недостатков:

• сложность теоретического исследования построенной модели и доказательства её сходимости к исходной задаче;
• Довольно невысокая точность и большая вычислительная сложность по сравнению с аналитическими методами решения данного класса задач.

 

Список использованных источников

1. А.А. Самарский Теория разностных схем. Учебное пособие. – М.: Наука, 1983. – 616 с.
2. Е.Н. Сметанникова Конспект лекций по курсу «Численные методы математического анализа». – Самара: Самарский государственный аэрокосмический университет. 1999. – 97с.
3. Дегтярев А.А. Метод конечных разностей. Электронное учебное пособие. – Самара: Самарский государственный аэрокосмический университет. 2011. – 83с.
4. Дегтярев А.А. Численные методы решения задач математической физики. – Самара: Самарский государственный аэрокосмический университет. 2009. – 26с.

 

Приложение

22 public class GridComputer {

23

24     private static int big_i = 10000;

25     private static int big_k = 10000;

26     private static double k = 0.067;

27     private static double c = 1.84;

28     private static double length = 3;

29     private static double big_t = 30;

 30     private static double x;

31     private static double t;

32

 

41     private static double maxDifference(double[][] numerical, int big_i, int big_k, double h_x, double t) {

42         double result = Double.NEGATIVE_INFINITY;

43         double current = 0;

44         SeriesComputer computer =

new SeriesComputer(length, k, c, 100500, big_i + 1, big_k + 1, 0, 0, big_t);

45         for (int i = 0; i < numerical[1].length; i++) {

46             current = Math.abs(numerical[1][i] - computer.computeInPoint(h_x * i, t, 30, big_t));

47             if (current > result){

48                 result = current;

 49             }

50         }

51         return result;

52     }

53

54     // Максимальное расхождение  сеточного и аналитического решения  для заданных параметров разбиения

 55     private static double eps(int big_i, int big_k) {

56         double result = Double.NEGATIVE_INFINITY;

57         double current;

58

59

60         double[][] currentLayer = new double[2][big_i + 1];

61         double h_t = big_t / big_k;

62         double h_x = length / big_i;

63         for (int i = 0; i < big_i + 1; i++) {

64             currentLayer[0][i] = i * h_x;

65             currentLayer[1][i] = initialFunction(i * h_x);

 66         }

67         double gamma = k / c * h_t / (h_x * h_x) * 0.5;

68         for (int k = 1; k <= big_k; k++) {

69             currentLayer[1] = computeOneLayer(currentLayer[1], gamma);

70             if (k > 1) {

71                 current = maxDifference(currentLayer, big_i, big_k, h_x, h_t * k);

72

73                 if (result < current) {

74                     result = current;

 75                 }

76             }

77         }

78

79

80         return result;

81     }

82

 

180

181     // значение t должно  принадлежать сетке

182     public static double[][] computeSeriesWhenT() {

183         double h_t = big_t / big_k;

184         long layer = Math.round(t / h_t);

185         double[][] result = computeLayer(layer);

186         return result;

187     }

188

189     public static double[][] computeSeriesWhenX() {

190         double h_t = big_t / big_k;

191         double h_x = length / big_i;

192         double[][] result = new double[2][big_k + 1];

193         int target_i = (int) Math.round(x / h_x);

194         double[][] currentLayer = new double[2][big_i + 1];

195         for (int i = 0; i < big_i + 1; i++) {

196             currentLayer[0][i] = i * h_x;

197             currentLayer[1][i] = initialFunction(i * h_x);

198         }

199         result[0][0] = 0;

200         result[1][0] = currentLayer[1][target_i];

201         double gamma = k / c * h_t / (h_x * h_x) * 0.5;

202         for (int k = 1; k <= big_k; k++) {

203             currentLayer[1] = computeOneLayer(currentLayer[1], gamma);

204             result[0][k] = k * h_t;

205             result[1][k] = currentLayer[1][target_i];

206         }

207         return result;

208     }

209

210     public static double initialFunction(double x) {

211         return Math.sin(Math.PI * x / length) + 10;

212     }

213

214     /**

215      * изображает трехдиагональную матрицу тремя массивами - диагоналями внутри

216      */

217     public static class ThreeDiagMatrix {

218

219         private double[] bot;

220         private double[] mid;

221         private double[] top;

222

223         public int size() {

224             return mid.length;

225         }

226

227         public double get(int i, int j) {

228             if (Math.abs(i - j) > 1) {

229                 return 0;

230             } else if (i == j) {

231                 return mid[j];

232             } else { // |i-j|=1

233                 if (i - j == 1) {

234                     return bot[j];

235                 } else { // (i - j == -1)

236                     return top[i];

237                 }

238             }

239         }

240

241         public void set(int i, int j, double value) {

242             if (Math.abs(i - j) > 1) {

243                 return;

244             } else if (i == j) {

245                 mid[j] = value;

246             } else { // |i-j|=1

247                 if (i - j == 1) {

248                     bot[j] = value;

249                 } else { // (i - j == -1)

250                     top[i] = value;

251                 }

252             }

253         }

254

255         public ThreeDiagMatrix(double[] bot, double[] mid, double[] top) {

256             if (!((bot.length == top.length) && (bot.length + 1 == mid.length))) {

257                 throw new java.lang.ArrayIndexOutOfBoundsException("Нельзя сделать трехдиагональную матрицу из этих массивов");

258             }

259             this.bot = Arrays.copyOf(bot, bot.length);

260             this.mid = Arrays.copyOf(mid, mid.length);

261             this.top = Arrays.copyOf(top, top.length);

262         }

263     }

264

265     /*

266      * решает систему  с трехдиагональной матрицей и столбцом справа b

267      */

268     public static double[] linsolve(ThreeDiagMatrix matrix, double[] b) {

269         int N = matrix.size() - 1;

270         // прямой  ход

271         // делим  первую строку

272         matrix.set(0, 1,

273                 matrix.get(0, 1) / matrix.get(0, 0));

274         b[0] = b[0] / matrix.get(0, 0);

275         matrix.set(0, 0, 1.0);

276         // делим остальные строки

277         for (int i = 1; i < matrix.size() - 1; i++) {

278             matrix.set(i, i + 1,

279                     matrix.get(i, i + 1) / (matrix.get(i, i) - matrix.get(i, i - 1) * matrix.get(i - 1, i)));

280             b[i] = (b[i] - matrix.get(i, i - 1) * b[i - 1]) / (matrix.get(i, i) - matrix.get(i, i - 1) * matrix.get(i - 1, i));

281             matrix.set(i, i - 1, 0);

282             matrix.set(i, i, 1.0);

283         }

284         // делим последнюю строку

285         b[N] = (b[N] - matrix.get(N, N - 1) * b[N - 1]) / (matrix.get(N, N) - matrix.get(N, N - 1) * matrix.get(N - 1, N));

286         matrix.set(N, N, 1.0);

287         matrix.set(N, N - 1, 0);

288         // обратный ход

289         double[] x = new double[matrix.size()];

290         x[N] = b[N];

291         for (int i = N - 1; i >= 0; i--) {

292             x[i] = b[i] - matrix.get(i, i + 1) * x[i + 1];

293         }

294         return x;

295     }

296

297     /*

298      * находит один  слой по времени с помощью  предыдущего. 

299      * gamma = a/2 * h_t / h_x^2 , где a это коэффициент при второй производной

300      */

301     public static double[] computeOneLayer(double[] preLayer, double gamma) {

302         int size = preLayer.length;

303         double[] top = new double[size - 3];

304         double[] mid = new double[size - 2];

305         double[] bot = new double[size - 3];

306         double[] b = new double[size-2];

307         for (int i = 0; i < size-3; i++) {

308             bot[i] = -1 * gamma;

309             mid[i] = (i==0) ? (1 + gamma) : (1 + 2 * gamma);

310             top[i] = -1 * gamma;

311             b[i] = (i == 0)

312                     ? ((1 - gamma) * preLayer[1] + gamma * preLayer[2])

313                     : ((1 - 2 * gamma) * preLayer[i+1] + gamma * (preLayer[i] + preLayer[i+2]));

314         }

315         mid[size - 3] = 1 + gamma;

316         b[size - 3] = ((1 - gamma) * preLayer[size-2] + gamma * preLayer[size-3]);

317         ThreeDiagMatrix matrix = new ThreeDiagMatrix(bot, mid, top);

318         double[] out = linsolve(matrix, b);

319         double[] result = new double[size];

320         for (int i=0; i<size-2; i++) {

321             result[i+1] = out[i];

322         }

323         result[0] = result[1];

324         result[size-1] = result[size-2];

325         return result;

326     }

327

328     public static double[][] computeLayer(long target_k) {

329         double[][] currentLayer = new double[2][big_i + 1];

330         double h_t = big_t / big_k;

331         double h_x = length / big_i;

332         for (int i = 0; i < big_i + 1; i++) {

333             currentLayer[0][i] = i * h_x;

334             currentLayer[1][i] = initialFunction(i * h_x);

335         }

336         double gamma = k / c * h_t / (h_x * h_x) * 0.5;

337         for (int k = 1; k <= target_k; k++) {

338             currentLayer[1] = computeOneLayer(currentLayer[1], gamma);

339         }

340         return currentLayer;

341     }

342

 

452 }