- •Часть I
- •Implicit none
- •Часть II Решение задачи на собственные значения для обыкновенного дифференциального уравнения второго порядка
- •Постановка задачи (вариант 2s).
- •Вывод разностной схемы для уравнения и граничных условий.
- •Данная система является 3-х диагональной,приведём подобные члены в уравнениях разностной схемы при -ых:
- •Implicit none
Минобрнауки России
Санкт-Петербургский государственный политехнический университет
Физико-технический факультет
Кафедра «Физика плазмы»
КУРСОВОЙ ПРОЕКТ
Численное решение краевых задач для обыкновенных дифференциальных уравнений второго порядка
Вариант - 2S
по дисциплине
«Численные методы и математическое моделирование»
Выполнил
студент гр.3101 В.А.Захаров
Руководитель
доцент, к.ф.-м.н. И.Ю.Веселова
«___» __________ 2012 г.
Санкт-Петербург
2012
Часть I
Решение краевой задачи для обыкновенного дифференциального уравнения второго порядка
Постановка задачи (вариант 2S).
Вывод разностной схемы для уравнения и граничных условий
Введём равномерную сетку :
Введём обозначение:
Используем приближённую формулу для вычисления интеграла:
И с помощью формулы центральных разностей получаем:
В итоге получим разностную схему для краевой задачи для обыкновенного дифференциального уравнения второго порядка.
Система из n+1 алгебраических уравнений:
Данная система является 3-х диагональной,приведём подобные члены в уравнениях разностной схемы при -ых:
Элементы нижней диагонали матрицы А - ,
главной диагонали - ,
верхней диагонали - ,
правой части -
Вывод выражения для главного члена погрешности аппроксимации уравнения и граничных условий.
Итак, мы свели краевую задачу к алгебраической задаче вида: Av=B, где A-матрица коэффициентов разностной схемы, В-вектор коэффициентов правой части, v-вектор точных решений системы Av=B
-вектор приближённых решений системы.
Общая погрешность решения краевой задачи складывается из погрешности аппроксимации и погрешности решения системы алгебраических уравнений. Рассмотрим подробнее погрешность аппроксимации уравнения.
, -невязка
Вычислим невязку :
Подставляем точное решение вместо приближённого и ищем невязку:
,
(1)
(2)
(3)
Тестовые задачи:
А)Тестовая задача с нулевой погрешностью аппроксимации:
Пусть u(r)=const=1 , k(r)=1 , q(r)=1
Подставим u,k,q в (1),(2),(3)
Тогда,получим что f(r)=1
Б)Тестовая задача с ненулевой погрешностью аппроксимации:
Пусть
Тогда используя (1),(2) и (3) получим,что
Код программы на фортране:
program Lab1
Implicit none
include 'link_fnl_static.h'
integer i
integer(8),parameter :: n=159,n1=n+1
real(8) k,q,f,u
external k,q,f,u
real(8) r(n1),delta,eps
real(8) c(n1),d(n1),e(n1),b(n1),nu,h
!nu=-60.0*exp(-4.0)
nu=0.0
h=2.0/n
r(1)=0.0
do i=2,n1
r(i)=r(i-1)+h
end do
open(1,file='out.dat')
write(1,*) "Количество узлов:",n1
write(1,*)
write(1,*) "Шаг разбиения:",h
write(1,*)
write(1,*) "Узлы основной сетки:"
do i=1,n1
write(1,*) r(i)
end do
write(1,*)
c(1)=0
e(n1)=0
!При i=1(апроксимация гр. условия 1 рода,r=0)
d(1)=-0.25*h*k(0.5*h)-q(r(1))*(h**3)/24
e(1)=0.25*h*k(0.5*h)
b(1)=-f(r(1))*(h**3)/24
!i=2,..,n (апроксимация во внутр. точках промежутка)
do i=2,n
c(i)=((r(i)-0.5*h)**2)*k(r(i)-0.5*h)/h
e(i)=((r(i)+0.5*h)**2)*k(r(i)+0.5*h)/h
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))
b(i)=-(r(i)**2)*f(r(i))*h
end do
!i=n+1 (апроксимация гр условия 2 рода,r=2)
c(n1)=((r(n1)-0.5*h)**2)*k(r(n1)-0.5*h)/h
d(n1)=-0.5*h*q(r(n1))*(r(n1)**2)-k(r(n1)-0.5*h)*((r(n1)-0.5*h)**2)/h
b(n1)=-0.5*h*f(r(n1))*(r(n1)**2)-nu*(r(n1)**2)
write(1,*) "Матрица коэффициентов:"
do i=1,n1
write(1,*) c(i),d(i),e(i),b(i)
end do
call DLSLTR(n1,c,d,e,b)
write(1,*)
write(1,*) "Приближенные и точные значения функции в узлах:"
write(1,*) "b(i): u(i):"
do i=1,n1
write(1,*)b(i),u(r(i))
end do
write(1,*)
write(1,*) "Погрешность апроксимации :"
write(1,*) "r(i): b(i)-u(r(i)):"
do i=1,n1
write(1,*) r(i),(b(i)-u(r(i)))
end do
delta=0.0 !Чебышевская норма погрешности
do i=1,n1
eps=abs(b(i)-u(r(i)))
if (eps>delta) then
delta=eps
end if
end do
write(1,*) "Чебышевская норма погрешности: ",delta
close(1)
end program Lab1
real(8) function k(r)
real(8) r
!k=r**2+1.0 !выданное задание
!k=r
k=1.0
return
end
real(8) function q(r)
real(8) r
q=1.0
return
end
real(8) function f(r)
real(8) r
!f=(-4*r**6+14*r**4+5*r**2-6)*exp(-r**2) !выданное задание
!f=60/exp(4.0)*r-7.5/exp(4.0)*r**2
f=1.0
return
end
real(8) function u(r)
real(8) r
!u=-7.5/exp(4.0)*r**2
u=1.0
return
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.