Минобрнауки России

Санкт-Петербургский государственный политехнический университет

Физико-технический факультет

Кафедра «Физика плазмы»

КУРСОВОЙ ПРОЕКТ

Численное решение краевых задач для обыкновенных дифференциальных уравнений второго порядка

Вариант - 2S

по дисциплине

«Численные методы и математическое моделирование»

Выполнил

студент гр.3101 В.А.Захаров

Руководитель

доцент, к.ф.-м.н. И.Ю.Веселова

«___» __________ 2012 г.

Санкт-Петербург

2012

Часть I

Решение краевой задачи для обыкновенного дифференциального уравнения второго порядка

1. Постановка задачи (вариант 2S).

2. Вывод разностной схемы для уравнения и граничных условий

Введём равномерную сетку :

Введём обозначение:

Используем приближённую формулу для вычисления интеграла:

И с помощью формулы центральных разностей получаем:

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

Система из n+1 алгебраических уравнений:

4. Данная система является 3-х диагональной,приведём подобные члены в уравнениях разностной схемы при -ых:

Элементы нижней диагонали матрицы А - ,

главной диагонали - ,

верхней диагонали - ,

правой части -

5. Вывод выражения для главного члена погрешности аппроксимации уравнения и граничных условий.

Итак, мы свели краевую задачу к алгебраической задаче вида: Av=B, где A-матрица коэффициентов разностной схемы, В-вектор коэффициентов правой части, v-вектор точных решений системы Av=B

-вектор приближённых решений системы.

Общая погрешность решения краевой задачи складывается из погрешности аппроксимации и погрешности решения системы алгебраических уравнений. Рассмотрим подробнее погрешность аппроксимации уравнения.

, -невязка

Вычислим невязку :

Подставляем точное решение вместо приближённого и ищем невязку:

,

(1)

(2)

(3)

6. Тестовые задачи:

А)Тестовая задача с нулевой погрешностью аппроксимации:

Пусть u(r)=const=1 , k(r)=1 , q(r)=1

Подставим u,k,q в (1),(2),(3)

Тогда,получим что f(r)=1

Б)Тестовая задача с ненулевой погрешностью аппроксимации:

Пусть

Тогда используя (1),(2) и (3) получим,что

7. Код программы на фортране:

1. program Lab1

2. Implicit none

3. include 'link_fnl_static.h'

4. integer i

5. integer(8),parameter :: n=159,n1=n+1

6. real(8) k,q,f,u

7. external k,q,f,u

8. real(8) r(n1),delta,eps

9. real(8) c(n1),d(n1),e(n1),b(n1),nu,h

10. !nu=-60.0*exp(-4.0)

11. nu=0.0

12. h=2.0/n

13. r(1)=0.0

14. do i=2,n1

15. r(i)=r(i-1)+h

16. end do

17. open(1,file='out.dat')

18. write(1,*) "Количество узлов:",n1

19. write(1,*)

20. write(1,*) "Шаг разбиения:",h

21. write(1,*)

22. write(1,*) "Узлы основной сетки:"

23. do i=1,n1

24. write(1,*) r(i)

25. end do

26. write(1,*)

27. c(1)=0

28. e(n1)=0

29. !При i=1(апроксимация гр. условия 1 рода,r=0)

30. d(1)=-0.25*h*k(0.5*h)-q(r(1))*(h**3)/24

31. e(1)=0.25*h*k(0.5*h)

32. b(1)=-f(r(1))*(h**3)/24

33. !i=2,..,n (апроксимация во внутр. точках промежутка)

34. do i=2,n

35. c(i)=((r(i)-0.5*h)**2)*k(r(i)-0.5*h)/h

36. e(i)=((r(i)+0.5*h)**2)*k(r(i)+0.5*h)/h

37. d(i)=-(r(i)+0.5*h)**2*k(r(i)+0.5*h)/h-(r(i)-0.5*h)**2*k(r(i)-0.5*h)/h-r(i)**2*h*q(r(i))

38. b(i)=-(r(i)**2)*f(r(i))*h

39. end do

40. !i=n+1 (апроксимация гр условия 2 рода,r=2)

41. c(n1)=((r(n1)-0.5*h)**2)*k(r(n1)-0.5*h)/h

42. d(n1)=-0.5*h*q(r(n1))*(r(n1)**2)-k(r(n1)-0.5*h)*((r(n1)-0.5*h)**2)/h

43. b(n1)=-0.5*h*f(r(n1))*(r(n1)**2)-nu*(r(n1)**2)

44. write(1,*) "Матрица коэффициентов:"

45. do i=1,n1

46. write(1,*) c(i),d(i),e(i),b(i)

47. end do

48. call DLSLTR(n1,c,d,e,b)

49. write(1,*)

50. write(1,*) "Приближенные и точные значения функции в узлах:"

51. write(1,*) "b(i): u(i):"

52. do i=1,n1

53. write(1,*)b(i),u(r(i))

54. end do

55. write(1,*)

56. write(1,*) "Погрешность апроксимации :"

57. write(1,*) "r(i): b(i)-u(r(i)):"

58. do i=1,n1

59. write(1,*) r(i),(b(i)-u(r(i)))

60. end do

61. delta=0.0 !Чебышевская норма погрешности

62. do i=1,n1

63. eps=abs(b(i)-u(r(i)))

64. if (eps>delta) then

65. delta=eps

66. end if

67. end do

68. write(1,*) "Чебышевская норма погрешности: ",delta

69. close(1)

70. end program Lab1

71. real(8) function k(r)

72. real(8) r

73. !k=r**2+1.0 !выданное задание

74. !k=r

75. k=1.0

76. return

77. end

78. real(8) function q(r)

79. real(8) r

80. q=1.0

81. return

82. end

83. real(8) function f(r)

84. real(8) r

85. !f=(-4*r**6+14*r**4+5*r**2-6)*exp(-r**2) !выданное задание

86. !f=60/exp(4.0)*r-7.5/exp(4.0)*r**2

87. f=1.0

88. return

89. end

90. real(8) function u(r)

91. real(8) r

92. !u=-7.5/exp(4.0)*r**2

93. u=1.0

94. return

95. end

8. Таблица результатов работы программы на тестах (зависимость Чебышёвской нормы погрешности решения от удваивающегося количества разбиений ,N=10).

k
1	1.110223024625157E-015	0,2
2	5.551115123125783E-015	0,335571
3	1.654232306691483E-014	0,207521
4	7.971401316808624E-014	0,488103
5	1.633138069223605E-013	----

Табл. 1(а).

k
1	2.518657457653639E-002	4,424779
2	5.700885613245221E-003	4,219409
3	1.355741532057197E-003	4,115226
4	3.306438260791689E-004	4,065041
5	8.147905095567953E-005	----

Табл. 1(б).

9. Оценка порядка точности разностной схемы найденного численного решения по табл. 1(а,б).

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

Видим, что с увеличением числа разбиений погрешность решения системы увеличивается. Это объясняется тем, что обусловленность матрицы растёт при росте количества узлов.

Из таблицы 1(б) мы видим, что погрешность уменьшается примерно в 4 раза при уменьшении шага в 2 раза, т.е. разностная схема имеет 2 порядок аппроксимации по h.

10. Таблица результатов работы программы для задачи, полученной от преподавателя.

	-0.02570	-0.00356	-0.00009	0,02214	0,00347	6,380403
	0.01373	0.03313	0.03722	0,0194	0,00409	4,743276
	0.10855	0.12977	0.13474	0,02122	0,00497	4,269618
	0.22533	0.24492	0.24963	0,01959	0,00471	4,159236
	0.31515	0.33203	0.33612	0,01688	0,00409	4,127139
	0.34802	0.36301	0.36667	0,01499	0,00366	4,095628
	0.32172	0.33639	0.33999	0,01467	0,0036	4,075
	0.25543	0.27099	0.27482	0,01556	0,00383	4,062663
	0.17548	0.19237	0.19653	0,01689	0,00416	4,060096
	0.10294	0.12098	0.12542	0,01804	0,00444	4,063063
	0.04829	0.06710	0.07173	0,01881	0,00463	4,062635

Табл. 2.

11. Оценка порядка точности разностной схемы найденного численного решения по табл. 2.

Из таблицы 2 видим, что разность численного решения в итых узлах при шаге h и h/2 примерно в 4 раза больше разности численного решения в тех же узлах при шаге h/2 и h/4, следовательно разностная схема имеет 2 порядок аппроксимации по h.

12. Выводы.

Практический результат согласуется с теорией, составленная разностная схема имеет второй порядок точности по h.