Лабораторная №7 Решение линейной краевой задачи методом прогонки

Лабораторная работа №7

Задание:

На отрезке [ a ; b ] решить методом прогонки линейную краевую задачу

y’’ + p(x)y’ +q(x)y = f(x)

C₁y(a) + C₂y’(a) = C

D₁y(b) + D₂y’(b) = D

α1 = 1,5 [ a ; b ] = [ 0 ; 0,6 ]

α2 = 2,5 C₁ = 1 C₂ = 0 C = 0,2

β = 33 D₁ = 1 D₂ = 0 D = 0,8

Псевдокод программы:

глобальные переменные скаляры N = 5

массив -- y[ N ] -- вещественный

прототипы подпрограмм – progonka ( float a[ N ] , float b[ N ] , float c[ N ] , float fi[ N ] )

алгоритм метод прогонки

начало

скаляры – I -- целое

-- a , b , c , d , c1 , c2 , d1 , d2 , h -- вещественные

массивы --al[ N ] , be[ N ] , ca[ N ] , fi[ N ] , x[ N ]

a = 0 b = 0.6

c = 0.2 d = 0.8

c1 = 1 d1 = 1

c2 = 0 d2 = 0

h = ( b – a ) / ( N – 1 )

be[ 0 ] = c1 * h - c2

ca[ 0 ] = c2

fi[ 0 ] = h * c

al[ N – 1 ] = - d2

be[ N – 1 ] = d1 * h + d2

fi[ N – 1 ] = d * h

для I = 0 ; I < N ; I = I+1

повторить

x[ I ] = a + I * h

для I = 0 ; I < N ; I = I+1

повторить

al[ I ] = 1 – h / 2 * ( ( 1.5 + 2.5 * x[ I ] ) / ( x[ I ] * x[ I ] – 1 ) )

be[ I ] = h * h * 33 / sqrt( 1 - x[ I ] * x[ I ] )

ca[ I ] = 1 + h / 2 * ( ( 1.5 + 2.5 * x[ I ] ) / ( x [ I ] * x[ I ] – 1 ) )

fi[ I ] = 0

вызов подпрограммы progonka ( al , be , ca , fi )

для I = 0 ; I < N ; I = I+1

повторить

вывод ( I , x[ I ] , I , y[ I ] )

Псевдокод подпрограммы:

алгоритм progonka ( float a[ N ] , float b[ N ] , float c[ N ] , float fi[ N ] )

начало

скаляры – I -- целое

-- e --вещественное

массивы -- u[ N ] , v[ N ]

u[ 0 ] = fi[ 0 ] / b[ 0 ]

v[ 0 ] = c[ 0 ] / b[ 0 ]

I = 1

для I = 0 ; I < N - 1 ; I = I+1

повторить

e = b[ I ] - a[ I ] * v[ I – 1 ]

u[ I ] = ( fi[ I ] - a[ I ] * u[ I – 1 ] ) / e

v[ I ] = c[ I ] / e

u[ N – 1 ] = ( fi[ N – 1 ] - a[ N – 1 ] * u[ N – 2 ] ) / ( b[ N – 1 ] - a[ N – 1 ] * v[ N – 2 ] )

y[ N – 1 ] = u [ N – 1 ]

для I = N – 2 ; I ≥ 0 ; I = I – 1

повторить

y[ I ] = u[ I ] - v[ I ] * y[ I + 1 ]

Запись программы на языке С++:

# include <stdio.h>

# include <conio.h>

# include <math.h>

# define N 5

float y[N];

progonka(float a[N],float b[N],float c[N],float fi[N]);

void main()

{ int i;

float a,b,c,d,c1,c2,d1,d2,h,al[N],be[N],ca[N],fi[N],x[N];

a=0; b=0.6; c=0.2; d=0.8;

c1=1; d1=1; c2=0; d2=0;

clrscr();

h=(b-a)/(N-1);

be[0]=c1*h-c2;

ca[0]=c2;

fi[0]=h*c;

al[N-1]=-d2;

be[N-1]=d1*h+d2;

fi[N-1]=d*h;

for(i=0;i<N;i++)

x[i]=a+i*h;

for(i=1;i<N-1;i++)

{al[i]=1-h/2*((1.5+2.5*x[i])/(x[i]*x[i]-1));

be[i]=h*h*33/sqrt(1-x[i]*x[i]);

ca[i]=1+h/2*((1.5+2.5*x[i])/(x[i]*x[i]-1));

fi[i]=0;}

progonka(al,be,ca,fi);

for(i=0;i<N;i++)

printf("\nx[%d]=%f y[%d]=%f",i,x[i],i,y[i]);

getch();

}

progonka(float a[N],float b[N],float c[N],float fi[N])

{ int i;

float e,u[N],v[N];

u[0]=fi[0]/b[0];

v[0]=c[0]/b[0];

i=1;

for(i=1;i<N-1;i++)

{

e=b[i]-a[i]*v[i-1];

u[i]=(fi[i]-a[i]*u[i-1])/e;

v[i]=c[i]/e;

}

u[N-1]=(fi[N-1]-a[N-1]*u[N-2])/(b[N-1]-a[N-1]*v[N-2]);

y[N-1]=u[N-1];

for(i=N-2;i>=0;i--)

y[i]=u[i]-v[i]*y[i+1];

}

Проверка программы ( тест ) :

Графики из полученных результатов:

При n = 5

x[0]=0.000000 y[0]=0.200000

x[1]=0.150000 y[1]=0.299735

x[2]=0.300000 y[2]=-0.530139

x[3]=0.450000 y[3]=0.070365

x[4]=0.600000 y[4]=0.800000

При n = 10

x[0]=0.000000 y[0]=0.200000

x[1]=0.066667 y[1]=0.661655

x[2]=0.133333 y[2]=-0.326648

x[3]=0.200000 y[3]=-0.697903

x[4]=0.266667 y[4]=0.487667

x[5]=0.333333 y[5]=0.735109

x[6]=0.400000 y[6]=-0.706514

x[7]=0.466667 y[7]=-0.771515

x[8]=0.533333 y[8]=1.032016

x[9]=0.600000 y[9]=0.800000

При n = 20

x[0]=0.000000 y[0]=0.200000

x[1]=0.031579 y[1]=-0.322512

x[2]=0.063158 y[2]=-0.199347

x[3]=0.094737 y[3]=0.346673

x[4]=0.126316 y[4]=0.198908

x[5]=0.157895 y[5]=-0.374276

x[6]=0.189474 y[6]=-0.198632

x[7]=0.221053 y[7]=0.406160

x[8]=0.252632 y[8]=0.198457

x[9]=0.284211 y[9]=-0.443455

x[10]=0.315789 y[10]=-0.198301

x[11]=0.347368 y[11]=0.487723

x[12]=0.378947 y[12]=0.198040

x[13]=0.410526 y[13]=-0.541192

x[14]=0.442105 y[14]=-0.197484

x[15]=0.473684 y[15]=0.607152

x[16]=0.505263 y[16]=0.196315

x[17]=0.536842 y[17]=-0.690671

x[18]=0.568421 y[18]=-0.193981

x[19]=0.600000 y[19]=0.800000