1

Treść zadania

Napisać program realizujący metodę Eulera dla zagadnienia :

x

(IV )

+ (1 + t

2

)x

000

+ (1 + 2t)x

00

+ tx = cos(1 + t) − 1

x(0) = x

0

(0) = x

00

(0) = x

000

(0) = 0

na przedziale [0,4].
Napisać program na rysowanie wykresów rozwiązań przybliżonych.

2

Teoria potrzebna do rozwiązania

Równanie liniowe n-tego rzędu
Mamy dane funkcje:
p

1

, p

2

, ..., p

n

: (a, b)− > R oraz f : (a, b)− > R

Oznaczamy przez y funkcje zmiennej t ∈ (a,b).
Liniowe rówanie różniczkowe n-tego rzędu:

(1) y

(n)

(t) + p

1

(t)y

(n−1)

(t) + ... + p

n−1

(t)y

0

(t) + p

n

(t)y(t) = f (t)

Niech t

0

∈ (a, b). Równanie spełnia warunek początkowy:

y(t

0

) = x

0

, y

0

(t

0

) = x

1

, ..., y

(n−1)

(t

0

) = x

n−1

Aby rozwiązać równanie definiujemy macierz A :

A =











0

1

0

. . .

0

0

0

1

. . .

0

. . .

. . .

. . .

. . .

. . .

0

0

0

. . .

1

−p

n

(t) −p

n−1

(t) −p

n−2

(t) . . . −p

1

(t)











oraz wektor : g(t) = [0, ..., 0, f (t)]

T

Niech x

∗

= (x

0

, x

1

, ..., x

n−1

)

T

.

Rozpatrujemy układ liniowy z warunkiem początkowym:

(2) z

0

(t) = A(t)z(t) + g(t) , z(t

0

) = x

∗

1

Twierdzenie
Jeśli φ : (a, b)− > R jest rozwiązaniem naszego równania różniczkowego z zagadnieniem
początkowym, to funkcja:
u(t) = (φ(t), φ

0

(t), ..., φ

(n−1)

(T ))

T

, t ∈ (a, b)

jest rozwiązaniem zagadnienia (2).
Jeśli u = (u

1

, ..., u

n

)

T

: (a, b)− > R

n

jest rozwiązaniem problemu (2), to φ = u

1

jest roz-

wiązaniem zagadnienia (1) z warunkiem początkowym.
Metoda Eulera
Niech f : [a, b]xR

n

− > R

n

będzie dana oraz η ∈ R

n

. Rozważamy zagadnienie Cauchy’ego:

(3) z

0

(t) = f (t, z(t)), z(a) = η

Niech h > 0 będzie ustalone oraz:

t

1

= a + ih, i = 0, 1, ..., N ,

a + N h = b

Zbiór t

0

, t

1

, ..., t

N

nazwamy siatką na przedziale [a, b].

Metoda Eulera:
y

0

= η

h

y

m+1

= y

m

+ hf (t

m

, y

m

) dla m = 0, 1, ..., N − 1

gdzie η

h

∈ R

n

jest dane.

Do wyznaczenia elementu y

m+1

korzystamy z elementu o jednym wcześniejszym nume-

rze, czyli y

m

oraz z warunku początkowego.

3

Działanie Programu

Program na początku liczy wartość h. Wartość ta zależy od wcześniej podanej liczby na-
turalnej N , tzn.

h = (B − A)/N ,
gdzie A = 0, B = 4.

Następnie wylicza wartości t

0

, t

1

, ..., t

N

(t

i+1

= t

i

+ h) i umieszcza je w tablicy t[].

Teraz korzystając z teori dotyczącej równań liniowych n-tego rzędu,
okazuje się, że program ma do policzenia układ czterech równań.

2

Dzialając według metody Eulera, program liczy wszystkie wartości i zapisuje je do czterech
tablic (z1,z2,z3,z4), wszystkie o długości N .

Na końcu, wyświetlane są wyliczone elementy tablicy węzłów (x-ów) oraz tablicy wyli-
czonych wartości funkcji będącej przybliżonym rozwiązaniem problemu (y-greków).

4

Wejście - wyjście

Wejście jest realizowane za pomocą funkcji z biblioteki simpio.h: GetInteger która zapew-
nia wystarczające zabezpieczenie przed wprowadzeniem błędnych danych. Wyjście jest
realizowane za pomocą standardowej funkcjiprintf, która w zupełności wystarcza.

5

Najważniejsza część kodu

void licz_wartosci()
{

int i;
z1[0]=z2[0]=z3[0]=z4[0]=0;
t[0]=A;
for (i=1;i<=N;i++)
{

t[i]=t[i-1]+h;
z1[i]=z1[i-1]+(h*z2[i-1]);
z2[i]=z2[i-1]+(h*z3[i-1]);
z3[i]=z3[i-1]+(h*z4[i-1]);
z4[i]=z4[i-1]+(h*((-1*t[i-1]*z1[i-1])+((-1-(2*t[i-1]))*z3[i-1])+((-1-(t[i-1]*t[i-1]))*z4[i-1]) + cos(1+t[i-1])-1));

}

}

W pętli for są obliczne kolejne wartości zerowej, pierwszej, drugiej i trzeciej pochodnej w
węzłach.

3

6

Przykładowe uruchomienia

6.1

Ekran

PROGRAM REALIZUJE METODE EULERA DLA ZAGADNIENIA:
x’’’’+(1+t^2)x’’’+(1+2t)x’’+tx = cos(1+t)-1 na przedziale [0,4]
Podaj N (calkowite)10

t[0] =0.000000

z1[0]=0.000000

t[1] =0.400000

z1[1]=0.000000

t[2] =0.800000

z1[2]=0.000000

t[3] =1.200000

z1[3]=0.000000

t[4] =1.600000

z1[4]=-0.011768

t[5] =2.000000

z1[5]=-0.062861

t[6] =2.400000

z1[6]=-0.189280

t[7] =2.800000

z1[7]=-0.411160

t[8] =3.200000

z1[8]=-0.716400

t[9] =3.600000

z1[9]=-1.088458

t[10]=4.000000

z1[10]=-1.516884

6.2

Plik wykres.txt

0.000000 0.000000

0.400000 0.000000

0.800000 0.000000

1.200000 0.000000

1.600000 -0.011768

2.000000 -0.062861

2.400000 -0.189280

4

2.800000 -0.411160

3.200000 -0.716400

3.600000 -1.088458

4.000000 -1.516884

5