Диф.ур 2 го порядка Опции

[-Andrey-]

Привиет народ, мне нужна задача:"Решение краевой задачи для дифура 2го порядка"
Нужно ее запрограммить в Паскале... может у кого она есть или ссылка на нее, то пожалуйста дайте....

[Altair]

Пусть на отрезке [p,q] задано диф. уравн. 2 порядка.

r(x) d^2y/dx^2 + s(x)dy/dx + t(x)y = f(x)

и краевые условия
(d1 dy/dx +g1y) |(x=p) = _f1
(d2 dy/dx +g2y) |(x=q) = _f2

функции r(x), s(x), t(x), f(x) заданны.
d1,d2,g1,g2,_f1,_f2 - известные числа.
требуется определить y=y(x).

Для численного решения, введем на отрезке [p,q] равномерную сетку:
x0 = p, xn = q,
h= (p-q)/n, x(k+1) - xk = h, k = 0,1.... n

h-шаг сетки

Для апроксимации производных в диф. уравнении будем использовать разностные соотношения

(dy/dx)xk = (yk+1 - yk-1)/2h
(d^2y/dx^2)xk = (yk+1 - 2yk + yk-1) / (h^2)

Для апросимации производных в краевых условиях будем использовать соотнощшения:
(dy/dx)p = (y1-y0)/h
(dy/dx)q = (yn-yn-1)/h

Окончательно, получим соедующие разностыне уравнения:
c0y0 - b0y1 = f0
-akyk-1 + ckyk - bkyk+1 = fk
-anyn-1 + cnyn = fn

здесь
c0 = g1 + d1/h, b0 = - d1/h
ak = S(xk)/2h - r(xk)/h^2, ck = -2r(xk)/h^2 +t(xk)
bk = - r(xk)/h^2 - s(xk)/2h, k = 1,...n-1
an=d2/h
cn=g2+d2/h


Записав систему в матричном виде, мы получим
Ax=f
где A будет трехдиагональная матрица...
и если (вроде обязательно), для нее будет выполненн критерий регулярности Адомара (то есть будут все параметры строчного преобладания положительны) то ее можно факторизовать
LU разложением, поулчим
LUx=f
теперь Ux=y
Ly=f это простая систама, и
Ux=y тоже...

на практике для 3-ч диаг, матрицы существует очень быстрый метод, в отличии от метода Гаусса - метод прогонки


вот программа на фортране, при желании можно транслировать в паскаль.

Код

Subroutine  progon (y,altha, beta, n, p, q)
    dimension y(1), altha(1), beta(1)
    h=(q-p)/n
    !h - шаг сетки
    x=p
    call coef(x,p,q,h,a,b,c,f)
    ! опр. коэфф. в уравнении
    altha(1)=b/c
    beta(1) = f/c
    do 5 i=1,n
        x=x+h
        call coef (x,p,q,h,a,b,c,f)
        beta(I+1) = (f+a*beta(i))/(c-a*altha(i))
        if (i.eq.n) go to 5
        altha(i+1) = b/(c-a*altha(i))
    5 continue
    y(n+1) = beta(n+1)
    !теперь прогоночные коэфф. определены
    do 10 i=n,1,-1
        y(i)=altha(i)*y(i+1)+beta(i)
    10 continue
    !вычислили знач. искомой функции
    return
end

Subroutine coef (x,p,q,h,a,b,c,f)
    call gran(d1,g1,fi1,d2,g2,fi2)
    !gran определяет коэфф. в краевых условиях
    if (x.eq.p) then
        a=0.0
        b=-d1/h
        c=g1+d1/h
        f=fi1
        return
    end if
    if (x.gt.q-h/2) then
        a=d2/h
        b=0.0
        c=g2+d2/h
        f=fi2
        return
    end if
    call func(x,r,s,t,f)
    ! func для точки Х опр. коэфф. R,S,T,F в исходной диф. уравнении
    a=s/(2*h)-r/(h**2)
    b=-r/(h**2) - s/(2*h)
    c=-2*r/(h**2)+t
    return
end

Subroutine  gran(d1,g1,fi1,d2,g2,fi2)
    d1=0.0
    g1=1.
    fi1=0.0
    d2=0.0
    g2=1
    fi2=0.0
    return
end

Subroutine Func (x,r,s,t,f)
    r=1.0
    s=0.0
    t=-2.0
    f=-3.0*sin(x)
    return
end


program TEST
    
    ! .......

    call progon (y,altha, beta, n, p, q)        
end program test

в ней, P и Q - границы отрезка... N-число узловых точек.
массив Y опр. значения искомой функции в узловых точках, массивы alpha и beta опр. прогоночные коэф..

В программе gran и func задаются исходные данные
в примере заданно уравнение
d^2/dx^2 - 2y = -3sinx
с граничными условиями Дирихле.

уфф..