Методы решения обыкновенных дифференциальных уравнений

Файл RK5.bas

function f(x as Double, y as Double) as Double

  f= exp(x*y)

end function

rem RK5

  dim x0 as Double,xN as Double

  dim h as Double

  dim y as Double

  dim k1 as Double, k2 as Double, k3 as Double, k4 as Double, k5 as Double

  dim x as Double

  input "Enter x0: ", x0

  input "Enter xN: ", xN

  input "Enter step: ", h

  input "Enter y(0): ", y

  x = x0

  print "x              y"

  do while x < xN+h

    print using "#.#              ##.######"; x; y

    k1 = h*f(x,y)/3

    k2 = h*f(x + h/3, y + k1)/3

    k3 = h*f(x + h/3, y + k1/2 + k2/2)/3

    k4 = h*f(x + h/2, y + 3*k1/8 + 9*k3/8)/3

    k5 = h*f(x + h, y + 3*k1/2 - 9*k3/2 + 6*k4)/3

    y = y + (k1 + 4*k4 + k5)/2

    x = x + h

  loop

end

Результат работы программы:

Метод Адамса

Файл adams.bas

function f(x as double,y as double) as double

  f = exp(x*y)

end function

rem adams

  rem dim y(10) as Double

  dim x as Double, y as Double ,y0 as Double

  dim x0 as Double,xN as Double

  dim h as Double,k1 as Double,k2 as Double

  dim k3 as Double, k4 as Double,k5 as Double

  dim i as Integer

  dim n as Integer

  input "Enter x0: ", x0

  input "Enter xN: " , xN

  input "Enter step: ", h

  input "Enter y(x0): ", y0

  n = abs(xN-x0)/h

  dim yh(0 to n) as Double

  dim yh2(0 to 2*n) as Double

  x = x0 : y = y0: i = 0

  do while x <= xN

    yh(i) = y

    k1 = h*f(x,y)

    k2 = h*f(x + h/2, y + k1/2)

    k3 = h*f(x + h/2, y + k2/2)

    k4 = h*f(x + h, y + k3)

    y = y + (k1+2*k2+2*k3+k4)/6

    x = x + h

    i = i + 1

  loop

  i = 4

  x = x0 + i*h

  do while x<=xN

    yh(i) = yh(i-1) + h*(55*f(x,yh(i-1))-59*f(x-h,yh(i-2))+37*f(x-2*h,yh(i-3))-9*f(x-3*h,yh(i-4)))/24

    x = x0 + i*h

    i = i + 1

  loop

  print "Result: "

  print using "y(#.#)= ##.######"; xN; yh(n)

  y = y0: x = x0: h = h/2: i = 0

  do while x <= xN

    yh2(i) = y

    k1 = h*f(x,y)

    k2 = h*f(x + h/2, y + k1/2)

    k3 = h*f(x + h/2, y + k2/2)

    k4 = h*f(x + h, y + k3)

    y = y + (k1+2*k2+2*k3+k4)/6

    x = x + h

    i = i + 1

  loop

  i = 8

  x = x0 + i*h

  do while x<=xN

    yh2(i) = yh2(i-1) + h*(55*f(x,yh2(i-1))-59*f(x-h,yh2(i-2))+37*f(x-2*h,yh2(i-3))-9*f(x-3*h,yh2(i-4)))/24

    x = x0 + i*h

    i = i + 1

  loop

  print using "eps = #.######"; (yh(n)-yh2(2*n))/(2^4-1)

end

Результат работы программы:

где eps – погрешность вычисленного решения.

Метод прогонки

Файл progon.bas

function funA(x as Double) as Double

  funA = 0.4+2*x+x*x/2

end function

function funB(x as Double) as Double

  funB = x^0.7

end function

function funC(x as Double) as Double

  funC = -0.95*x*x

end function

function funF(x as Double) as Double

  funF = (sqr(x)+cos(2*x/3.1416))/funA(x)

end function

REM PROGON

  dim i as Integer

  dim a as Double, b as Double, x as Double, h as Double

  dim NS as Integer

  input "Enter N: ", NS

  dim y(0 to NS+1) as Double

  dim c(0 to NS+1) as Double

  dim d(0 to NS+1) as Double

  dim m(0 to NS+1) as Double

  dim n(0 to NS+1) as Double

  input "Enter a: ", a

  input "Enter b: ", b

  if a>=b then

    print "incorrect interval"

    end

  end if

  print using "Enter y(#.##)"; a

  input "", y(0)

  print using "Enter y(#.##)"; b

  input "", y(NS)

  h = (b-a)/NS

  m(0) = -2+funB(a)*h/funA(a)

  n(0) = 1-funB(a)*h/funA(a) + funC(a)*h*h/funA(a)

  c(0) = 1/m(0)

  d(0) = -n(0)*y(0)+funF(a)*h*h

  for i = 1 to NS

    x=a+i*h

    m(i) = -2+funB(x)*h/funA(x)

    n(i) = 1-funB(x)*h/funA(x)+funC(x)*h*h/funA(x)

    c(i) = 1/(m(i)-n(i)*c(i-1))

    d(i) = funF(x)*h*h-n(i)*c(i-1)*d(i-1)

  next i

  for i = NS to 2 step (-1)

    y(i-1) = c(i-2)*(d(i-2)-y(i))

  next i

  for i=0 to NS

    print using "y(#.###)= ##.######"; a+i*h; y(i)

  next i

end

Результат работы программы:

Сводные таблицы результатов

Таблица 1. Результаты решения ОДУ

Таблица 2. Результаты решения краевых задач

Вывод

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

Полученные результаты, как видно из сводной таблицы, имеют некоторые расхождения в силу особенностей языков программирования. Результаты решения ОДУ методами Рунге-Кутты 4 и 5 порядков, Адамса одинаковы, так как данные методы имеют четвертый порядок точности, в отличие от метода Эйлера, который имеет меньший порядок точности (первый порядок точности).

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

Список используемой литературы

1. Демидович, Марон «Основы вычислительной математики», М. , 1970.
2. Алексеев Е.Р., Чеснокова О.В. Пакет Scilab. Начало работы.
3. Конспект лекций по Фортрану – 2001.

39