Ministerul Educaţiei, Culturii şi
Cercetării al Republicii Moldova
Universitatea Tehnică a Moldovei
Departamentul Informatică și
Ingineria Sistemelor
RAPORT
despre lucrarea
de laborator nr.7
la Metode Numerice
Tema: Integrarea
numerică a ecuațiilor
diferențiale
Varianta 7
A îndeplinit: [Link]. CR-221fr, Ciobanu Stanislav
A controlat: [Link]., Seiciuc V.
Chişinău, 2023
Scopul lucrării
1) Să se determine soluţiile numerice ale ecuaţiilor diferenţiale pe segmentul [ a , b ] prin
metoda Euler, Euler modificat, și Runge-Kutta cu pasul h=0.05.
2) Să se efectueze o analiză a rezultatelor primite.
Realizarea lucrărilor
Conform variantei mele, cu numărul 7, am următoarea problemă (ecuații diferențiale)
dată spre rezolvare:
' cos y 2
y= −0.25 y , y ( 0 )=0; a=0 ; b=1 ;
1.25+ x
În prima parte a sarcinii, care este mai mult teoretică, se cere de descris metodele
utilizate, astfel problema Cauchy pentru ecuaţia diferenţială de ordinul n de forma
y (n )=f ( x , y , y ' , K , y n−1 ) constă în determinarea funcţiei y= y (x), care satisface ecuaţia dată și şi
condiţiile iniţiale y ( x 0 )= y 0 , y '(x ¿¿ 0)= y 01 , … , y ( x 0 )= y 0 ,n−1 ¿ , unde x 0 , y 0 , … , y 0 ,n−1 se
n−1
consideră numere cunoscute. Modelele matematice a multor procese în tehnică conţin sisteme
de ecuaţii diferenţiale ordinare de forma:
{
d y1
=f 1 (t , y 1 , y 2 , K , y n )
dt
d y2
=f 2 (t , y 1 , y 2 , K , y n )
dt
KKKKKKKKKKKKKK
d yn
=f n (t , y 1 , y 2 , K , y n)
dt
În sistemul dat funcţiile f i ( t , y 1 , y 2 ,… , y n ) (i=1 ,n) sunt cunoscute. Pentru sistemul dat,
problema Cauchy permite determinarea funcţiilor y i ( t ) (i=1, n) care satisfac toate ecuaţiile
sistemului şi condiţiile iniţiale y 1 ( t 0 ) = y 10, y 2 ( t 0 ) = y 20, ..., y n ( t 0 )= y n 0. În cazul în care se poate
construi soluţia generală a ecuaţiei diferențiale de ordinul n sau sistemului dat, problema se
reduce la determinarea constantelor de integrare, care satisfac relaţiile respective. În majoritatea
cazurilor determinarea soluţiei pentru problema Cauchy este imposibilă ceea ce impune construirea
soluţiei aproximative. Metodele aproximative în dependenţă de modul în care se obţine soluţia
problemei Cauchy pot fi împărţite în două grupe. Prima este cea a metodelor analitice, care oferă
soluţia sub forma unei expresii analitice. Cea de-a doua este a metodelor numerice care permit
obţinerea soluţiei aproximative sub forma unui tabel de valori ale soluţiei într-un şir de puncte.
Astfel, metoda seriilor de puteri (metoda derivării succesive) este o metodă analitică.
V-om considera problema Cauchy, presupunînd că pentru ecuaţie se îndeplinesc condiţiile de
existenţă şi unicitate a soluţiei. Fie că soluţia particulară y= y (x) a ecuaţiei permite dezvoltarea ei în
serie Taylor în jurul punctului x=x 0:
' '' (n)
y ( x0 ) y ( x0 ) y ( x0 )
y ( x )= y ( x 0 )+ ( x−x 0 ) + ( x−x 0 ) 2+ Λ + ( x−x 0 )n + Λ
1! 2! n!
(i )
Condiţiile iniţiale permit ca valorile y ( x 0 ) (i =0 , n−1) să fie cunoscute. Pe de altă parte
(n )
ecuaţia se poate determina y ( x 0 ) . Ceilalţi coeficienţi din dezvoltarea Taylor se vor calcula prin
diferenţierea succesivă a ecuaţiei, substituind în ultima valorile cunoscute deja a derivatelor de
ordin superior în punctul x=x 0. În mod similar se poate determina soluţia problemei Cauchy.
De exemplu, să se determine primii 5 termeni ai dezvoltării în serie a soluţiei y= y (x) pentru
2
ecuaţia y ' ' ( 1+ y )=( y ' ) + y ' cu condiţiile iniţiale y ( 0 )=2 și y ' ( 0 )=2. V-om căuta soluţia sub
forma seriei de puteri:
y' ( 0) y '' ( 0 ) 2 y (n ) ( 0 ) n
[ ]
'
'' y
y ( x )= y ( 0 ) + x+ x + Λ+ x + Λ y ( 0 )= (1+ y ') =2
1! 2! n! 1+ y x =0
Fie funcţia f (x , y ) satisface în dreptunghiul D=¿{|x−x 0|≤ a ;| y− y 0|≤ b ; condițiile
|f ( x1 , y 1 )−f ( x1 , y 2 )|≤ N |y 1− y 2|
df df
|
= +f
df
dx dx dy |
≤ M (N , M −const )
hM
Atunci eroarea metodei se estimează în modul următor | y ( x n )− y n|≤
n
[ ( 1+hN ) −1], unde
2N
y ( x n ) este soluția exactă, iar soluția y n numerică în punctul x=x n, obţinută la pasul n . Această
estimare este o estimare teoretică, care se utilizează rar. În practică se procedează în felul
următor: se determină soluţia după metoda Euler cu pasul h , iar apoi cu pasul h /2 . Atunci
eroarea poate fi estimată conform relaţiei | y ( x n )− y |≈| y n− y n |, ( y n este soluția numerică în
¿ ¿ ¿
punctul x=x n pentru pasul egal cu h /2).
Haideți să vorbim și despre metoda Euler modificată. După cum s-a văzut eroarea
metodei este considerabilă, care poate fi întrucîtva micşorată dacă se va folosi un pas de
discretizare destul de mic. Aceasta din urmă cere eforturi de calcul mari. De aceea au fost
propuse mai multe modificări ale metodei Euler în scopul ridicării ordinului preciziei. Una din
aceste modificări constă în precizarea iteraţională a soluţiei. Fie că este cunoscută soluţia y i
numerică a problemei anterioare. Considerând această soluţie ca o aproximaţie iniţială,
determinăm următoarea aproximaţie:
y i+ 1= y i+ hf (x i , y i )
( )
' '' '2
' ( 1+ y ) y − y
' ' ''
y y y (
+ 1+ y )
' '' ' '' '
y = ( y +1) = 2
; y ( 0 )=2
1+ y 1+ y ( 1+ y )
IV 1
y = ¿
( 1+ y )3
IV
y ( 0 )=2
Aşadar soluţia aproximativă se poate scrie sub două forme:
2 2 2 3 2 4
y ( x ) ≈ 2+2 x+ x + x + x
2! 3! 4!
2 1 3 1 4
y ( x ) ≈ 2+2 x+ x + x + x
3 12
Metoda Euler este o metodă numerică care permite determinarea soluţiei aproximative y (x )
a problemei Cauchy sub forma unui tabel de valori. Ideea metodei se bazează pe dezvoltarea funcţiei
y= y (x) în seria Taylor în vecinătăţile unui şir de puncte, numite noduri x=x i (i=0 , 1 ,2 , …) şi din
care se elimină termenii seriei ce conţin derivatele de ordin mai superior ca unu. Să considerăm
problema Cauchy pentru ecuaţia diferenţială de ordinul unu y ' =f ( x , y ) și y ( x 0 )= y 0 . Să
considerăm un pas de discretizare destul de mic h şi construim un şir de puncte (noduri)
echidistanţate x i=x 0 +ih(i=0 , 1, 2 , …). Soluţia aproximativă a problemei se obţine conform
algoritmului y i+ 1= y i+ hf ( xi , y i ) (i=0 , 1 , 2, …).
Geometric, metoda Euler constă în înlocuirea curbei y= y (x) prin linia poliginală ABC
(ca în desenul de mai jos) în modul următor: tangenta la curba integrală y=0 în punctul
A(x 0 , y 0 ) are coeficientul unghiular y 0=f ' (x 0 , y 0) . Punctele B şi C se obţin în rezultatul
soluţionării problemei Cauchy prin metoda Euler. După fiecare pas noi de fapt trecem pe o altă
curbă integrală. Astfel segmentul BC este de acum un segment al tangentei dusă prin punctul
B(x 1 , y 1 ) la curba integrală 1. Deci sensul geometric al metodei constă în înlocuirea curbei
integrale O care trece prin punctul A(x O , y O ) printr-o linie poligonală, numită frînta lui Euler.
Metoda lui Euler este o metodă directă, care foloseşte informaţia legată de o singur pas al
algoritmului, adică din acest punct de vedere este o metodă unipas. Eroarea metodei se
constituie din două componente, și anume din eroarea legată de trunchierea seriei Taylor și
eroarea rotungirii. Deoarece aceste erori se adună pas cu pas, eroarea soluţiei într-un punct x ,
aflat la o distanţă l de la punctul x O, poate fi destul de mare. Din acest punct de vedere, metoda
Euler este o metodă de ordinul 1 de exactitate. Folosind această valoare se obţine algoritmul
modificat Euler:
h
y i+ 1= y i+
2
[ f ( x i , y i) + f (x +1 , y i+1 )]
Metoda Runge-Kutta este de asemenea o metodă directă de rezolvare a ecuaţiilor
diferenţiale ordinare. Ideea metodei constă în construirea unei formule de calcul al valorilor
soluţiei problemei Cauchy în punctele x i (i=0 , 1 ,2 , …) de tipul y i+ 1= y i+ hφ(x i , y i , h) în care
funcţia φ (x i , y i , h) ar aproxima segmentul seriei Taylor în jurul punctului x i cu exactitatea 0 (
p +1
h ) şi totodată să nu conţină derivatele funcţiei f (x , y ). Trebuie de menţionat că metoda
Runge-Kutta de ordinul 1 ( p=1 ) este de fapt metoda lui Euler. Vom exemplifica pentru cazul
metodei Runge-Kutta de ordinul 2 cum poate fi obţinută formula de recurenţă. Vom căuta
funcţia φ sub forma φ ( x , y , h )=Af ( x , y ) + Bf (x+ αh , y + βh) în care A , B,α , β sunt nişte
constante temporar necunoscute. Pentru a le calcula vom dezvolta termenul al doilea din partea
dreaptă a expresiei în seria Taylor în împrejurimile punctului x i, păstrând numai primii trei
termeni (exactitatea de aproximare 0(h2)). Obţinem:
φ ( x , y , h )=Af ( x i , y i ) + [ ∂ f ( xi , yi )
∂x
hα +
∂ f ( xi , y i )
∂y
2
hβf ( x i , y i ) + 0(h ) ]
Pe de altă parte soluţia y (x ) a problemei, fiind dezvoltată în seria Taylor cu exactitatea 0(h2 )
este:
[ h ' x,y ' x,y
y i+ 1= y i+ h f ( x i , yi ) + (f x( )+ f y( ) f ( x i , y i ) ) +0 (h2)
2
i i i i
]
Comparând ultimele două formule obţinem un sistem de ecuaţii pentru calculul coeficienţilor
necunoscuţi A+ B=1; αB=1 /2 ; βB=1/2. Sistemul dat cu trei ecuaţii conţine patru necunoscute
şi este un sistem nedeterminat, soluţia căreia este A=1−λ ; B= λ ; α=β=γ /2 λ ;(λ ≠ 0) şi care
depinde de paramentrul λ . În particular pentru λ=1/2 se obţine formula Runge-Kutta de ordinul
2:
{
h
y i+1= y i + ( k 1 +k 2 ) , unde(i=0 ,1 , 2)
2
k 1=f ( x i , yi ) ; k 2=f ( x i +h , y i+ h k 1 )
Într-un mod similar poate fi obţinută schema de calcul al metodei Runge-Kutta de orice ordin
p. În particular metoda Runge-Kutta de ordinul patru constă în aplicarea algoritmului:
{
y i +1= y i + Δ y i
1 (i) (i) i=0 ,1 , 2
Δ y=
6
( k 1 +2 k 2 +2 k 3 +k 4 )
(i) (i)
unde:
(i) (i) (i)
k 1 =hf ( x i , y i ) ; k 2 =hf (xi + h/ 2; y i + k 1 /2)
(i) (i) (i) (i )
k 2 =hf ( xi +h /2; y i +k 2 / 2); k 3 =hf (x i+ h ; y i+ k 3 )
Pasul h poate fi schimbat de la punct la punct. Pentru a verifica corectitudinea alegerii pasului
h se poate calcula mărimea:
| |
(i) (i )
k 2 −k 3
θ=
k 1(i) −k 2(i)
care nu trebiue să depăşească cîteva sutimi, în caz contrar el trebuie micşorat. Metoda
Runge-Kutta are ordinul h 4 de exactitate pe întreg segmentul de calcul a soluţiei. Estimarea erorii
metodei este destul de dificilă, însă această estimare poate fi efectuată după regula Runge în modul
¿
următor. Fie y (x n ) valoarea soluţiei exacte în punctul x n, iar y n și y n sunt valorile aproximative ale
soluţiei în acest punct, calculate cu pasul h /2 şi h . Atunci eroarea metodei este determinată de relaţia:
¿
Concluzie
În cadrul acestei lucrări de laborator, care a avut un caracter deja mai mult practic, am
însușit trei metode (Euler, Euler modificat și Runge-Kutta) pentru calculul numeric ale ecuațiilor
diferențiale, obținând aproximativ aceeleași soluții indiferent de formula aleasă. Metoda Euler,
Euler modificat și metoda Runge-Kutta sunt toate metode numerice utilizate pentru rezolvarea
ecuațiilor diferențiale ordinare (EDO).
Metoda Euler este mai puțin precisă și poate conduce la erori mari, în special pentru pași
mari. Pe când metoda Euler modificată oferă o precizie mai bună decât metoda Euler simplă.
Metoda Runge-Kutta oferă o precizie superioară, de obicei, mai precisă decât metodele Euler.
Din punct de vedere al complexității computaționale, metoda Runge-Kutta este mai complexă,
dar oferă o bună balanță între precizie și eficiență. Metoda Euler (simplă) poate fi nesigură în
anumite situații, în special pentru ecuații cu condiții inițiale critice. Metodele Euler modificat și
Runge-Kutta sunt mai stabile și pot gestiona o varietate mai mare de condiții inițiale fără să
dezvolte instabilități. În general, alegerea între aceste metode depinde de specificul problemei pe
care o rezolvați și de echilibrul pe care doriți să îl atingeți între precizie și eficiență
computațională. Metoda Runge-Kutta de ordinul patru (RK4) este adesea preferată datorită
echilibrului său între precizie și complexitate.
Vreau să menționez că deși am auzit de la colegi că au utilizat Wolfram Alpha pentru a se
verifica, eu totuși am decis să caut calculatoare specializate online pentru calculul numeric al
integralelor. Astfel am găsit platforma [Link], care prin intermediul plugin-ului javascript
MathJax, m-a ajutat să calculez soluțiile specifice fiecărei din cele 3 metode, ca ulterior să le
compar cu rezultatele obținute de programul meu. Din punct de vedere al programului vreau doar
să spun că nu m-am ciocnit cu dificultăți, la fel am creat o interfață plăcută cu meniu de
interacționare om-mașină.
Cursul de metode numerice face parte din disciplinele fundamentale de pregătire a
studenţilor din domeniul ingineriei, având ca scop prezentarea principiilor şi relaţiilor de calcul
matematic numeric care stau astăzi la baza costrucţiei programelor de calcul profesionale
utilizate în prezent de orice inginer. De aceea, cu siguranță voi căuta domenii de aplicare a
acestuia în sfera mea de activitate profesională.
Anexe
Anexa 1. Listing-ul programului C++
(cu screenshot-urile rezultatelor)
// Includem bibliotecile necesare
#include <iostream>
#include <cmath>
using namespace std;
// Definirea ecuatiei diferentiale
double f(double x, double y) {
return cos(y)/(1.25+x)-0.25*pow(y,2);
// Metoda Euler
void EULER(double a, double b, double h) {
cout<<"x\t\t y\n";
cout<<"----------\t ----------\n";
double x=a;
double y=0.0;
while(x<=b) {
cout<<x<<"\t\t "<<y<<"\n";
y+=h*f(x,y);
x+=h;
// Metoda Euler modificat
void EULER_MOD(double a, double b, double h) {
cout<<"x\t\t y\n";
cout<<"----------\t ----------\n";
double x=a;
double y=0.0;
while(x<=b) {
cout<<x<<"\t\t "<<y<<"\n";
double k1=h*f(x,y);
double k2=h*f(x+h,y+k1);
y+=0.5*(k1+k2);
x+=h;
// Metoda Runge-Kutta (de ordinul 4)
void RUNGE_KUTTA(double a, double b, double h) {
cout<<"x\t\t y\n";
cout<<"----------\t ----------\n";
double x=a;
double y=0.0;
while(x<=b) {
cout<<x<<"\t\t "<<y<<"\n";
double k1=h*f(x,y);
double k2=h*f(x+0.5*h,y+0.5*k1);
double k3=h*f(x+0.5*h,y+0.5*k2);
double k4=h*f(x+h,y+k3);
y+=(k1+2*k2+2*k3+k4)/6;
x+=h;
// Functia programului principal
int main() {
system("clear");
// Declaram limitile intervalului
double a=0.0;
double b=1.0;
// Declaram valoarea pasului
double h=0.05;
int opt;
cout<<"Ecuatia diferentiala data spre rezolvare pe intervalul [a="<<a<<",
b="<<b<<"]"<<endl;
cout<<"cu pasul h="<<h<<" este y'=cos(y)/(1.25+x)-0.25*y^2\n"<<endl;
cout<<"------------------------------------------------------------------"<<endl
;
cout<<"| Meniul cu lista metodelor de rezolvare:
|"<<endl;
cout<<"------------------------------------------------------------------\
n"<<endl;
cout<<"1. Metoda numerica Euler\n";
cout<<"2. Metoda Euler modificata\n";
cout<<"3. Metoda directa Runge-Kutta\n";
cout<<"0. Iesire din program\n"<<endl;
cout<<"Selectati optiunea: ";
cin>>opt;
switch(opt){
case 1:
system("clear");
cout<<"--------------------------------------------------------"<<endl;
cout<<"| Rezultatele metodei numerice Euler de rezolvare |"<<endl;
cout<<"--------------------------------------------------------\
n"<<endl;
EULER(a,b,h);
cout<<endl;
getchar();
exit(1);
break;
case 2:
system("clear");
cout<<"--------------------------------------------------------"<<endl;
cout<<"| Rezultatele metodei Euler modificate |"<<endl;
cout<<"--------------------------------------------------------\
n"<<endl;
EULER_MOD(a,b,h);
cout<<endl;
getchar();
exit(1);
break;
case 3:
system("clear");
cout<<"--------------------------------------------------------"<<endl;
cout<<"| Rezultatele metodei directe Runge-Kutta |"<<endl;
cout<<"--------------------------------------------------------\
n"<<endl;
RUNGE_KUTTA(a,b,h);
cout<<endl;
getchar();
exit(1);
break;
default:
cout<<"Eroare! Ai ales o optiune gresita!";
break;
case 0:
return 0;