Metode directe Runge-Kutta. Pentru simplitate, luați în considerare un spațiu bidimensional al variabilelor x și y și o mulțime deschisă G care îi aparține. Fie definită pe această mulțime deschisă o funcție diferențiabilă continuu f(x, y) și dată o ecuație. L

Rezolvarea problemei Cauchy pentru sistemul ODE prin metoda Runge-Kutta de ordinul 4
Algoritm secvenţial
Complexitate secvenţială 4mn
Cantitatea de date de intrare m + 3
Volumul datelor de ieșire (m+1)n
Algoritm paralel
Înălțimea unei forme paralele în etaje Pe)
Lățimea formei paralele cu etaje O(m)

Autorii principali ai descrierii: A.A. Forester (1.1,1.2,1.3,1.4,1.5,1.6,1.7,1.8,1.9), D.A. Alimov (1.1,1.7,1.10,2.4,2.7)

1 Proprietăți și structura algoritmului

1.1 Descrierea generală a algoritmului

Metoda Runge-Kutta de ordinul al patrulea este cea mai comună metodă din familia metodei Runge-Kutta de algoritmi numerici pentru rezolvarea ecuațiilor diferențiale obișnuite și a sistemelor acestora. Aceste metode iterative pentru aproximări explicite și implicite au fost dezvoltate în jurul anului 1900 de către matematicienii germani K. Runge și M. V. Kutta.

Formal, metoda Runge-Kutta este o metodă Euler modificată și corectată, sunt scheme de ordinul doi de precizie. Există scheme standard de ordinul al treilea, care nu sunt utilizate pe scară largă.
Ideea principală a algoritmilor Runge-Kutta este de a înlocui partea dreaptă a ecuației diferențiale, care depinde de funcția necunoscută dorită, cu o anumită aproximare. Dacă problema Cauchy este rescrisă în formă integrală

y(x) = y_(0) + \int\limits_(x_(0))^(x)f(t,y)dt,

apoi, prin aplicarea diverselor formule numerice pentru calcularea integralei din partea dreaptă a ecuației, se pot obține metode Runge-Kutta de diverse ordine.

Vedere generală a formulelor metodelor Runge - Kutta cu pasul de grilă h_(n) :

y = f(t,y),\,\,\, y(t_(0)) = y_(0), y_(n+1) = y_(n) + h_(n+1)\sum\limits_(i=1)^(s)b_(i)K_(n)^(i), K_(n)^(i) = f(t_(n) + c_(i)h_(n+1), y_(n) + h_(n+1)\sum\limits_(j=1)^(i -1)a_(ij)K_(n)^(j)).

Metoda specifică Runge-Kutta este determinată de un set de coeficienți b_(i), c_(j), a_(ij), care trebuie să satisfacă anumite relaţii.

1.2 Descrierea matematică a algoritmului

1.2.1 Metoda Runge-Kutta de ordinul 4 pentru problema Cauchy pentru DE de ordinul I

Luați în considerare problema Cauchy, unde partea dreaptă satisface condiţiile teoremelor de existenţă şi unicitate pentru soluţie.

y" = f(x,y),\ a \leq x \leq b;\ y(a) = y^0

Configurați o grilă uniformă

x_i = a + ih,\ i = 1,\dots, n,\ h = \frac(b-a)(n)

Să introducem notația y(x_i) = y_i . Obtinem formula de calcul:

\begin(cases) k_1 = hf(x_i,y_i)\\ k_2 = hf(x_i + h/2,y_i + k_1/2)\\ k_3 = hf(x_i + h/2,y_i + k_2/2)\ \ k_4 = hf(x_i + h,y_i + k_3)\\ y_(i+1) = y_i + [ k_1 + 2k_2 + 2k_3 + k_4 ]/6 \\ \end(cases)

1.2.2 Metoda Runge-Kutta de ordinul 4 pentru problema Cauchy pentru un sistem DE de ordinul întâi

Rezolvarea numerică a problemei Cauchy pentru sistemele de EDO de ordinul I prin metodele Runge-Kutta se caută prin aceleași formule ca și pentru EDO de ordinul I. De exemplu, o soluție Runge-Kutta de ordinul 4 poate fi găsită punând:

y_i \rightarrow \bar y_i f(x_i,y_i) \rightarrow \bar f(x_i,\bar y_i) k_l \rightarrow \bar k_l \bar k_l = \begin(pmatrix) k^i_(l,1)\\ \vdots\\ k^i_(l,m)\\ \end(pmatrix)

unde m este dimensiunea sistemului, l = 1, \dots, 4. Drept urmare, obținem

\begin(cases) \bar k_1 = h\bar f(x_i,\bar y_i)\\ \bar k_2 = h\bar f(x_i + h/2,\bar y_i + \bar k_1/2)\\ \ bara k_3 = h\bar f(x_i + h/2,\bar y_i + \bar k_2/2)\\ \bar k_4 = h\bar f(x_i + h,\bar y_i + \bar k_3)\\ \ bar y_(i+1) = \bar y_i + [ \bar k_1 + 2\bar k_2 + 2\bar k_3 + \bar k_4 ]/6 \\ \end(cases)

1.3 Miezul de calcul al algoritmului

1.4 Macrostructura algoritmului

Macrostructura algoritmului constă în calcularea coeficienților k_(j),\,j=1,..,4 la aflarea valorilor funcția dorită la fiecare nod.

1.5 Schema de implementare a algoritmului secvenţial

Iată o posibilă implementare a algoritmului secvenţial în limbajul Matlab (partea dreaptă a ODE este dată ca funcţie func(x,y))

1 funcție = RungeKutta4 (y0, a, b, n) 2 h = (b - a ) / n ; 3 x = a : h : b 4 5 y (:, 1 ) = y0 ; 6 pentru i = 2 : n 7 k1 = h * func (x (i ), y (:, i )) 8 k2 = h * func (x (i ) + h / 2 , y (:, i ) + k1 / 2 ) 9 k3 = h * func (x (i ) + h / 2 , y (:, i ) + k2 / 2 ) 10 k4 = h * func (x (i ) + h , y (:, i ) + k3 ) 11 y (:, i + 1 ) = y (:, i ) + (k1 + 2 * k2 + 2 * k3 + k4 ) / 6 12 sfârşitul 13 sfârşitul

1.6 Complexitatea secvenţială a algoritmului

Fiecare pas al buclei algoritmului constă din 4 apeluri la funcția \bar f, 11 înmulțiri și 10 adunări. Deoarece apelul la funcția \bar f este cea mai complexă operație, complexitatea execuției liniare a algoritmului poate fi definită ca 4mn .

1.7 Graficul de informare

Să descriem graficul algoritmului analitic și grafic.
Să fie dat un sistem de m ecuații diferențiale unidimensionale. Graficul are o structură tridimensională, care este planuri (nivele) mărginite interconectate. La fiecare i-lea nivel, se calculează i-a aproximare a vectorului soluție al sistemului. Există n astfel de niveluri în total, unde n este numărul de noduri în care se calculează soluția aproximativă a sistemului de ecuații diferențiale. Ieșirea nivelului i este intrarea la nivelul i+1.
Vârfurile graficului de informații sunt împărțite în două tipuri:

  • Primul tip de vârfuri corespunde calculului unei coordonate a funcției care se află în partea dreaptă a ecuației diferențiale, în punctele furnizate la intrarea vârfului. Rezultatul operației este o coordonată vectorială k_(l),\,\,\,l = 1,2,3,4. Vor exista 4m de astfel de vârfuri la fiecare nivel, iar numărul lor total este de 4mn.
  • Operația a + bc corespunde celui de-al doilea tip de vârfuri. Aceste vârfuri pot fi, de asemenea, împărțite în două grupuri. La vârfurile primului grup se calculează expresii, al căror rezultat va fi dat ca argumente pentru vârfurile primului tip, în timp ce datele de intrare pentru aceste vârfuri vor fi rezultatele declanșării a m vârfuri de primul tip situate pe acelasi nivel. Vor fi 3m de vârfuri la fiecare nivel al primului grup, iar numărul lor total va fi de 3mn. Datele de intrare pentru vârfurile celui de-al doilea grup vor fi rezultatul declanșării unui vârf de primul tip și unul de același vârf de al doilea tip. Rezultatul declanșării acestui vârf la nivelul i-lea va fi m coordonate ale vectorului i-lea aproximare a soluției sistemului. Vor exista 4m de astfel de vârfuri la fiecare nivel, iar numărul lor total este de 4mn.

Să oferim o ilustrare grafică a graficului de informații (Fig. 1). Pentru a nu aglomera prea mult graficul, luați în considerare cazul unui sistem de două ecuații diferențiale. Vârfurile primului tip sunt marcate cu roșu, vârfurile celui de-al doilea tip sunt marcate cu verde. Datele de intrare ale problemei sunt transmise la toate vârfurile primului tip, precum și la m vârfuri din stânga ale celui de-al doilea grup al celui de-al doilea tip și la toate nodurile primului grup al celui de-al doilea tip de la primul nivel, sunt furnizate valorile inițiale ale funcției dorite (vezi paragraful ). Datele de ieșire vor fi rezultatele declanșării a m vârfuri drepte ale celui de-al doilea grup de al doilea tip la fiecare nivel. Rezultatele declanșării vârfurilor rămase sunt datele intermediare ale algoritmului.

Figura 1. Graficul de informații pentru un sistem de două ecuații diferențiale.

1.8 Resursa de paralelism al algoritmului

Deoarece în cele de mai sus circuit de calcul cea mai laborioasă este operația de calcul a părților corecte ale ODE la calcul k_i (i = 1, \dots, 4), atunci atenția principală trebuie acordată paralelizării acestei operații. Aici se va aplica abordarea descompunerii ecuațiilor sistemului ODE în subsisteme. Prin urmare, pentru inițializare, luați în considerare următoarea schemă de descompunere a datelor pentru elementele procesorului disponibile cu memorie locală: pentru fiecare \mu - PE (element de procesor) ( \mu = 0, \dots, p-1) este distribuit n/p ecuații diferențiale și vectorul \bar y_0 . Calculele suplimentare se fac conform următoarei scheme:

  1. pe fiecare PE, n/p componente corespunzătoare ale vectorului \bar k_1 sunt calculate simultan prin formula [ \bar k_1 ]_(\mu) = h[ \bar f(x_i, \bar y_i) ]_(\mu)
  2. pentru a asigura cea de-a doua etapă de calcul este necesară asamblarea integrală a vectorului \bar k_1 pe fiecare PE. Apoi, componentele vectorului \bar k_2 sunt calculate independent prin formula [ \bar k_2 ]_(\mu) = h[ \bar f(x_i + h/2,\bar y_i + \bar k_1/2)]_(\mu);
  3. pe fiecare PE se realizează ansamblul vectorului \bar k_2, se calculează componentele vectorului \bar k_3:\ [ \bar k_3 ]_(\mu) = h [\bar f(x_i + h/2,\bar y_i + \bar k_2/2)]_(\mu);
  4. pe fiecare PE se realizează ansamblul vectorului \bar k_3, se calculează componentele vectorului \bar k_4:\ [ \bar k_4 ]_(\mu) = h [\bar f(x_i + h,\bar y_i + \bar k_3)]_(\mu);
  5. componentele vectoriale sunt calculate cu paralelism ideal \bar y_(i+1):\ [\bar y_(i+1)]_(\mu) = [\bar y_(i)]_(\mu) + ([ \bar k_1 ]_(\mu) ) + 2[ \bar k_2 ]_(\mu) + 2[ \bar k_3 ]_(\mu) + [ \bar k_4 ]_(\mu))/6\ iar vectorul \bar y_(i+1) este asamblat pe fiecare PE. Dacă trebuie să continuați proces de calcul, apoi setăm i = i + 1 și trecem la articolul 1

Rețineți că acest algoritm are paralelism finit, dar nu paralelism de masă, deoarece ciclurile algoritmului sunt dependente de informații.

La fiecare nivel acest algoritm fiecare PE efectuează patru operații pentru calcularea părților corecte ale ODE, șaisprezece operațiuni pentru adăugarea de vectori și înmulțirea unui vector cu un număr, precum și transferul de date p-1 între alte PE, ceea ce încetinește destul de mult algoritmul și este peste cap. la paralelizarea algoritmului. În acest caz, numărul de iterații buclei este egal cu lungimea vectorului x, adică n. Din cele de mai sus rezultă că la clasificarea după înălțimea JPF, algoritmul are complexitatea O(n) , iar la clasificarea după lățimea JPF, este O(m) (cu condiția ca p = m ).

1.9 Date de intrare și de ieșire ale algoritmului

Datele de intrare ale algoritmului sunt:

  1. vector y^0 de dimensiunea m ;
  2. limitele intervalului de timp a și b ;
  3. rata de eșantionare n ;

Dimensiunea totală de intrare: m + 3

Ieșirea este

  1. n vectori \bar y_i de dimensiunea m ;
  2. vector \bar x de dimensiunea n

Dimensiunea totală a ieșirii: (m+1)n

1.10 Proprietăţile algoritmului

Enumerăm principalele proprietăți ale algoritmului:

2 Implementarea software a algoritmului

2.1 Caracteristici ale implementării algoritmului secvenţial

2.2 Localitatea datelor și calculului

2.2.1 Implementarea locală a algoritmului

2.2.1.1 Structura acceselor la memorie și evaluarea calitativă a localității
2.2.1.2 Cuantificare localitate
2.2.1.3 Analiza Apex-Map

2.3 Metode și caracteristici posibile ale implementării paralele a algoritmului

2.4 Scalabilitate a algoritmului și implementarea acestuia

2.4.1 Scalabilitate a algoritmului

2.4.2 Scalabilitatea implementării algoritmului

Să studiem scalabilitatea implementării paralele a metodei Runge-Kutta conform metodologiei. Studiul a fost realizat pe supercomputerul Lomonosov al Complexului de Supercomputer al Universității din Moscova. Am studiat propria noastră implementare paralelă a algoritmului, scrisă în C++ folosind standardul MPI. Programul a fost compilat folosind un compilator de Intel versiunea 15.0 (fără a specifica cheile de compilare) folosind biblioteca OpenMPI versiunea 1.8.4. Calculele au fost efectuate pe noduri din grupul T-Blade2 cu procesoare Intel Xeon 5570 Nehalem 2,93 GHz, cu un nucleu alocat pentru fiecare proces MPI.

Setul și limitele valorilor parametrilor variabili pentru lansarea implementării algoritmului:

  • număr de procesoare în trepte de puteri de la două până la 64, apoi în trepte de 10;
  • dimensiunea sistemului.

În urma experimentelor, s-a obținut următorul interval de eficiență a implementării algoritmului:

  • eficiență minimă de implementare 0,0000032%;
  • eficiența maximă de implementare este de 0,0024%.

Iată câteva caracteristici ale implementării paralele testate:

  • pentru testare s-a folosit o parte dreaptă pseudo-aleatorie a sistemului, constând din suprapuneri de funcții elementare selectate aleatoriu;
  • valorile inițiale ale funcției dorite au fost stabilite aleatoriu;
  • întrucât în ​​cazul general partea dreaptă a sistemului este necunoscută, nu este posibilă măsurarea performanței în gigaflopi; numărul de puncte la care a fost evaluată funcția a fost folosit ca măsură de performanță.
  • Figura 3. Eficiența algoritmului paralel în funcție de dimensiunea sistemului și de numărul de procesoare.

    Să construim estimări de scalabilitate pentru implementarea paralelă testată a metodei Runge-Kutta de ordinul 4:

    • După numărul de procese -0,000000768. Odată cu creșterea numărului de procese, eficiența în intervalul considerat de modificări ale parametrilor de lansare scade, dar scăderea este destul de lentă, ceea ce se datorează eficienței maxime extrem de scăzute a implementării paralele a algoritmului. Scăderea eficienței cu creșterea numărului de procese pentru sistemele cu o dimensiune mică se explică prin faptul că, dacă numărul de procese depășește dimensiunea sistemului, unele dintre ele nu vor participa la calcule. Pentru sistemele cu dimensiuni mari, creșterea proceselor afectează negativ eficiența, deoarece numărul de transferuri între ele crește, ceea ce încetinește semnificativ munca.
    • După dimensiunea sarcinii 0,000000763. Odată cu creșterea dimensiunii sistemului cu un număr fix de procese, crește eficiența în intervalul considerat de modificări ale parametrilor de lansare, iar creșterea se observă în aproape toată regiunea luată în considerare. Acest lucru se explică prin faptul că, cu un număr fix de procese, o creștere a dimensiunii sistemului duce la o creștere a sarcinii de calcul pe fiecare procesor, în urma căreia timpul crește. munca de calcul referitor la timpul de nefuncţionare. Odată cu creșterea numărului de procese, rata de creștere a eficienței cu creșterea dimensiunii problemei scade, deoarece numărul de transferuri devine suficient de mare.
    • În două direcții 0,00000000174. Odată cu creșterea simultană a numărului de procese și a dimensiunii sistemului, eficiența crește. Rata de creștere a eficienței este extrem de lentă, ceea ce se explică prin costuri generale ridicate.

    2.5 Performanța dinamică și eficiența implementării algoritmului

    2.6 Concluzii pentru clasele de arhitectură

    2.7 Implementări existente ale algoritmului

    Metoda de ordinul 4 este cea mai des folosită dintre toate schemele Runge-Kutta, așa că există multe implementări software în serie ale acesteia, atât comerciale, cât și gratuite. Una dintre cele mai cunoscute implementari software este ode45 in Mediul MATLAB, care are un numar mare de oportunități pentru configurarea unui calcul numeric, ceea ce îl face destul de convenabil de utilizat. De asemenea, metoda Runge-Kutta de ordinul 4 este implementată în biblioteca paralelă PETSc, care este distribuită liber. În ciuda faptului că în această bibliotecă majoritatea metodelor sunt implementate folosind algoritmi paraleli, metoda Runge-Kutta din ea este implementată secvenţial. Este destul de dificil de găsit o implementare paralelă a metodei în cazul luării în considerare sarcină comună Cauchy, dar sunt multe lucrări științifice, care descriu implementări paralele ale metodei Runge-Kutta pentru clase specifice de sisteme, de exemplu, în munca lui A.V. Starchenko descrie o implementare paralelă pentru sisteme liniare.

Metoda clasică Runge-Kutta 4 ordine

Subiectul 6.3. Metoda Runge-Kutga

Metode Runge-Kutta- o familie importantă de algoritmi numerici pentru rezolvarea (sistemelor) de ecuații diferențiale obișnuite. Aceste metode iterative de calcul aproximativ explicit și implicit au fost dezvoltate în jurul anului 1900 de către matematicienii germani K. Runge și M. V. Kutta.

Formal, metodele Runge-Kutta sunt o metodă Euler modificată și corectată, sunt scheme de ordinul doi de precizie. Există scheme standard de ordinul al treilea, care nu sunt utilizate pe scară largă. Cel mai des folosit și implementat în diverse pachete matematice (Maple, MathCAD, Maxima) schema standard al patrulea ordin (vezi mai jos). Uneori când se efectuează calcule cu precizie crescută sunt utilizate scheme de ordinul al cincilea și al șaselea.

Metoda Runge-Kutta de ordinul 4 este atât de răspândită încât este adesea menționată pur și simplu ca metoda Runge-Kutta.

Luați în considerare problema Cauchy. Apoi valoarea din următorul punct se calculează cu formula:

Valoarea pasului de grilă pentru .

Această metodă are ordinul 4, adică eroarea la fiecare pas este O(h5), iar eroarea totală la intervalul final de integrare este O(h4).

Familia metodelor directe Runge-Kutta este o generalizare a metodei Runge-Kutta de ordinul 4. Este dat de formule

Metoda specifică este determinată de numărul s și de coeficienții bi,aij și ci. Acești coeficienți sunt adesea aranjați într-un tabel

Să fie dată o ecuație diferențială de ordinul întâi cu condiția inițială y(x 0)=y 0. Alegem pasul h și introducem notația:

x i = x 0 + ih și y i = y(x i), unde i = 0, 1, 2, ... .

Soluția este similară cu metoda descrisă mai sus.

ecuație diferențială. Diferența constă în împărțirea pasului în 4 părți.

Conform metodei Runge-Kutta de ordinul al patrulea, valorile succesive y i ale funcției dorite y sunt determinate de formula:

yi+1 = yi+?yi

unde i = 0, 1, 2...

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

numerele k1, k2 ,k3, k4 la fiecare pas sunt calculate prin formulele:

k1 = h*f(x i, y i)

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

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

k4 = F(x i +h, y i +k3)*h

Aceasta este o metodă explicită în patru pași cu 4 ordine de precizie.

Diagrama bloc a procedurii de rezolvare a unei ecuații diferențiale prin metoda Runge-Kutta este prezentată în Figura 6.

F(x, y) - o funcție dată - trebuie descrisă separat.

Parametrii de intrare:

X0, XK - inițial și final

independent

variabil;

Y0 - valoarea y 0 din condiția inițială

N este numărul de segmente de partiție;

Parametri de ieșire:

Y - matrice de valori ale soluției dorite

la nodurile grilei.

Algoritm pentru rezolvarea problemei

Soluție în MathCAD

Lista de programe este activată Limbajul vizual De bază

Dim xr(), yr(), xe(), ye(), xo(), yo() Ca Single

Privat x0, y0, h, xk, k1, k2, c, k3, k4, yd Ca singur

Privat n, i Ca întreg

Funcția publică f(ByVal a, ByVal b) ca unică

f = -(a - b) / a

Private Sub Button1_Click(ByVal sender As System.Object, ByVal e As System.EventArgs) Se ocupă de Button1.Click

x0 = Val(TextBox1.Text)

xk = Val(TextBox2.Text)

y0 = Val(TextBox4.Text)

h = Val(TextBox3.Text)

n = (xk - x0) / h

c = y0 / x0 + Math.Log(x0)

DataGridView1.ColumnCount = 4

DataGridView1.RowCount = n + 2

DataGridView1.Item(0, 0).Value = "x"!}

DataGridView1.Item(1, 0).Value = "General"!}

DataGridView1.Item(2, 0).Value = "Euler M"!}

DataGridView1.Item(3, 0).Value = "Runge Kutt"!}

Pentru i = 0 până la n - 1

xe(i) = Math.Round((xe(0) + i * h), 2)

tu(i + 1) = tu(i) + h * f(xe(i) + h / 2, tu(i) + h / 2 * f(xe(i), tu(i)))

DataGridView1.Item(2, 1).Value = ye(0)

DataGridView1.Item(0, 1).Value = xe(0)

DataGridView1.Item(0, i + 1).Valoare = xe(i)

DataGridView1.Item(2, i + 1).Value = Str(ye(i))

Pentru i = 0 până la n - 1

xr(i) = Math.Round((xe(0) + i * h), 2)

k1 = h * f(xr(i), yr(i))

k2 = h * f(xr(i) + h / 2, yr(i) + k1 / 2)

k3 = h * f(xr(i) + h / 2, yr(i) + k2 / 2)

k4 = h * f(xr(i) + h, yr(i) + k3)

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

yr(i + 1) = yr(i) + yd

DataGridView1.Item(3, 1).Value = yr(0)

DataGridView1.Item(3, i + 1).Value = Str(yr(i))

Pentru i = 0 până la n - 1

xo(i) = Math.Round((xe(0) + i * h), 2)

yo(i) = xo(i) * (c - Math.Log(xo(i)))

DataGridView1.Item(1, 1).Value = yo(0)

DataGridView1.Item(1, i + 1).Value = Str(yo(i))

Chart1.Series.Add(„Soluție generală”)

Chart1.Series(„Soluție generală”).ChartType = SeriesChartType.Line

Pentru i = 0 până la n - 1

Diagrama 1.Seria(„Soluție generală”).Puncte.AddXY(xo(i), yo(i))

Chart1.Series(„Soluție generală”).ChartArea = „ChartArea1”

Chart1.Series.Add(„Euler M”)

Chart1.Series("Euler M").ChartType = SeriesChartType.Point

Pentru i = 0 până la n - 1

Diagrama 1.Seria(„Euler M”).Puncte.AddXY(xe(i), ye(i))

Chart1.Series("Euler M").ChartArea = "ChartArea1"

Chart1.Series("Euler M").Culoare = Culoare.Albastru

Chart1.Series.Add(„Runge Kutt”)

Chart1.Series(„Runge Kutt”).ChartType = SeriesChartType.Line

Pentru i = 0 până la n - 1

Diagrama 1.Seria(„Runge Kutt”).Puncte.AddXY(xr(i), yr(i))

Chart1.Series("Runge Kutt").ChartArea = "ChartArea1"

Chart1.Series(„Runge Kutt”).Culoare = Culoare.Verde

si informatica

Institutul Tehnic de Comunicații și Informatică din Ural

Departamentul de Fizică, Matematică Aplicată și Informatică.

LUCRARE DE CURS

in informatica:

Vizualizarea metodelor numerice.

Rezolvarea ecuațiilor diferențiale obișnuite.

Completat de: Arapov

Gr. ME-71

Verificat de: Minina E.E.

Ekaterinburg

2008

Introducere………………………………………………………………….3

1. Enunțul problemei…………………………………………………….4

2. Descrierea metodelor de soluţionare……………………………………………………..5

2. 1. Esența sarcinii……………………………………………………………………….5

2. 2. Sensul geometric al problemei……………………………………………….5

2. 3. Metode numerice de rezolvare a problemei Cauchy…………….6

2. 4. Metoda Euler ……………………………..…………………………….6

2. 5. Metoda Runge-Kutta de ordinul al 4-lea……………………………………………….8

2. 6. Rezolvarea problemei enunţate de către Euler şi

Runge-Kutta de ordinul al 4-lea ……………………………………………….9

2. 6. 1. Metoda Euler ………………….. ……….……………9

2. 6. 2. Metoda Runge-Kutta de ordinul al 4-lea ……………………………10

3. Algoritmul de rezolvare a problemei…………………………………………...11

3. 1. Algoritmi de subrutine.………………………………………....11

3. 1. 1. Subrutina metodei Euler…………..………..11

3. 1. 2. Subrutina metodei Runge-Kutta de ordinul al 4-lea ………..12

3. 1. 3. Subrutină soluție generală ……….. …………………1 3

1 3

3. 3. Algoritmul programului…………………………………………………………………14

4. Forma programului…………………………………………………….17

5. Lista programelor…………………………………………………………..18

6. Rezolvarea problemei în MathCad …………………………………………..20

Concluzie………………………………………………………………………22


Introducere

Multe legi fizice, care sunt supuse anumitor fenomene, sunt scrise sub forma unei ecuații matematice care exprimă o anumită relație între unele mărimi.Importanța mare pe care o au ecuațiile diferențiale pentru matematică, și mai ales pentru aplicațiile ei, se explică prin faptul că studiul multor materiale fizice și sarcini tehnice. Ecuațiile diferențiale joacă un rol esențial în alte științe, cum ar fi biologia, economia și ingineria electrică; de fapt, ele apar oriunde este nevoie de o descriere cantitativă (numerică) a fenomenelor. Prin urmare, rezolvarea ecuațiilor diferențiale va fi întotdeauna o sarcină necesară și relevantă.

Acest termen de hârtie este soluția unei ecuații diferențiale prin două metode numerice: metoda Euler și metoda Runge-Kutta a 4 ordine de precizie.

Pentru a atinge acest obiectiv, mi-am propus următoarele sarcini:

  1. Scrieți un program pentru a rezolva ecuația diferențială dată prin două metode numerice din program Visual Basic.
  2. Verificați soluția cu aplicația MathCad.
  3. Comparați primit metode diferite rezultate din moment ce solutie comuna.


1. Enunțarea problemei

Rezolvați problema Cauchy pentru o ecuație diferențială de ordinul I pe segmentul [ X 0; X k ] cu pasul h si starea initiala: Y (X 0 ) = Y 0 .

Răspunsul ar trebui primit sub forma unui tabel cu rezultate:

Y (1)

Y (2)

Y 0 (1)

Y 0 (2)

Y(X0)

Y 1 (1)

Y 1 (2)

Y(X1)

Y k (1)

Y k (2)

Y(X k )

Unde Y (1) , Y (2) soluții obținute prin diverse metode numerice, YT rezolvarea exactă a unei ecuații diferențiale.

Este posibil să se prezinte rezultatele soluției nu sub formă de tabel, ci sub formă de liste.

Datele tabelului sunt vizualizate pe formular sub formă de grafice.

Înainte de a calcula ultima coloană a tabelului de rezultate, este necesar să se scadă valoarea coeficientului din condițiile inițiale c utilizat în soluția generală.

Ecuație diferențială

Decizie comună

y = - x y/(x+1)

1 ,2

0, 1

y=c (x+1) exp(-x)


2. Descrierea metodelor de soluţionare

2. 1. Esența problemei

Pentru a rezolva o ecuație diferențială obișnuită, este necesar să se cunoască valorile variabilei dependente și (sau) derivatele acesteia pentru unele valori ale variabilei independente. Dacă acestea termeni suplimentari sunt date pentru o valoare a variabilei independente, atunci o astfel de problemă se numește problemă cu condiții inițiale sau problemă Cauchy. Timpul este adesea folosit ca variabilă independentă în problema Cauchy.

Problema Cauchy poate fi formulată astfel:

Să fie date ecuația diferențială și condiția inițială y (x 0) = y 0 . Trebuie să găsim funcția y( X ) care satisface atât ecuația indicată, cât și condiția inițială.

Rezolvarea numerică a problemei Cauchy se reduce la tabelarea funcției dorite.

Graficul soluției unei ecuații diferențiale se numește curbă integrală.

2. 2. Sensul geometric al problemei

y = f (x, y ) este tangenta pantei tangentei la graficul soluției în punctul (x, y) la axa 0X, este panta (Fig. 1).

Figura 1. Sensul geometric al problemei Cauchy.

Existența soluției:

Dacă partea dreaptă f (x, y ) este continuă într-un anumit domeniu R , definit de inegalități

| x - x0 |< а; | y - y 0 | < b ,

atunci există cel puțin o soluție y = y(x) definită în vecinătatea |x x 0 | < h , где h este un număr pozitiv.

Această soluție este unică dacă R condiția Lipschitz este îndeplinită

| f (x , y )- f (x , y )| ≤ N | y - y |(x, y),

unde N este o constantă (constanta Lipschitz) în funcție, în cazul general, de a și b. Dacă f(x , y) are o derivată mărginită

f y (x, y) în R , atunci putem seta N = max | f y (x, y)| pentru (x, y ) aparținând lui R .

2. 3. Metode numerice de rezolvare a problemei Cauchy

Când se utilizează metode numerice, segmentul [х 0, X ] - zone de schimbare continuă a argumentului x de către mulțime. constând dintr-un număr finit de puncte x 0 < х 1 < ... < x n = Х - сеткой.

În acest caz x i sunt numite noduri de grilă.

Multe metode folosesc grile uniforme cu un pas:

Problema Cauchy definită mai devreme pe intervalul continuu [x 0, X ], este înlocuit cu analogul său discret - un sistem de ecuații, prin rezolvarea căruia se pot găsi succesiv valorile y 1 , y 2 ,…, y n - valorile aproximative ale funcției la nodurile grilei.

Rezolvarea numerică a problemei Cauchy este utilizată pe scară largă în diverse zoneștiință și tehnologie, iar numărul de metode dezvoltate pentru aceasta este destul de mare. Aceste metode pot fi împărțite în următoarele grupuri.

Metode într-un singur pas în care să găsiți următorul punct
curba y = f(x ) necesită informații doar despre un singur pas anterior.
Metodele cu un singur pas sunt metoda Euler și metodele Runge-Kutta.

Metode de prognoză și corecție (multi-step), în care să se găsească următorul punct al curbei y = f(x ) necesită informații despre mai mult de unul dintre punctele anterioare. Pentru a obține o valoare numerică suficient de precisă, se recurge adesea la iterație. Printre astfel de metode se numără metodele lui Milne, Adams-Bashfort și Hamming.

Metode explicite de care funcția Ф nu depinde y n +1.

Metode implicite de care depinde funcția Ф y n +1.

2.4 Metoda Euler.

Uneori, această metodă este numită metoda Runge-Kutta de ordinul întâi de precizie.

Aceasta metoda un pas. Tabularea funcției are loc alternativ în fiecare punct. Pentru a calcula valoarea funcției la următorul nod, trebuie să utilizați valoarea funcției la un nod anterior.

Să fie dată o ecuație diferențială de ordinul întâi:

Y = f(x, y)

cu starea initiala

y (x 0 ) = y 0

Să alegem un pas h și introduceți notația:

x i \u003d x 0 + ih și y i \u003d y (x i), unde i \u003d 0, 1, 2, ...,

x i - noduri de grilă,

y eu - valoarea functiei integrale in noduri.

Ilustrațiile pentru soluție sunt prezentate în Figura 2.

Desenați o dreaptă AB printr-un punct ( x i, y i ) la un unghi α. în care tan α = f (x i , y i )

În conformitate cu semnificația geometrică a problemei, linia AB este tangentă la funcția integrală. Să înlocuim punctul funcției integrale cu un punct situat pe tangentei AB.

Atunci y i +1 = y i + Δy

Echivalează părțile drepte tg α = f (x i , y i ) și. obține

Prin urmare Δ y \u003d h ∙ f (x i, y i).

Înlocuiți în această expresie formula yi +1 = yi + ∆y , apoi transformați-l. Ca rezultat, obținem formula de calcul al următorului punct al funcției integrale:

Figura 2. Metoda Euler.

Din formula se poate observa că pentru a calcula fiecare punct următor al funcției integrale este necesar să se cunoască valoarea unui singur punct anterior. Astfel, cunoscând condițiile inițiale, se poate construi o curbă integrală pe un interval dat.

Diagrama bloc a procedurii de rezolvare a unei ecuații diferențiale prin metoda Euler este prezentată în Figura 3.

F(x , y) - funcția dată trebuie

fi descris separat.

Parametrii de intrare:

X0, XK initial si final

valorile variabilei independente;

Y 0 valoare y 0 din starea initiala

y (x 0 ) = y 0 ;

Parametri de ieșire:

Y - matrice de valori ale soluției dorite

la nodurile grilei.

Figura 3. Diagrama bloc a procedurii de rezolvare a unei ecuații diferențiale folosind metoda Euler.

Metoda Euler este una dintre cele mai simple metode de rezolvare numerică a ecuațiilor diferențiale obișnuite. Dar dezavantajul său semnificativ este o eroare mare de calcul. În Figura 2, eroarea de calcul pentru i -r o pasul este notat cu ε. Cu fiecare pas, eroarea de calcul crește.

2.5 Metoda Runge-Kutta 4 ordine

Să fie dată o ecuație diferențială de ordinul întâi cu condiția inițială y (x 0 )= y 0. Alegeți pasul h și introduceți notația:

x i = x 0 + ih și y i = y (x i ), unde i = 0, 1, 2, ... .

Soluția este similară cu metoda descrisă mai sus.

ecuație diferențială. Diferența constă în împărțirea pasului în 4 părți.

Conform metodei Runge-Kutta de ordinul al patrulea, valori succesive y i a funcției dorite y sunt determinate de formula:

y i+1 = y i +∆y i unde i = 0, 1, 2 ...

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

a numere k 1, k 2 , k 3, k 4 la fiecare pas sunt calculate prin formulele:

k 1 \u003d h * f (x i, y i)

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

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

k4 = F(x i + h, y i + k3) * h

Aceasta este o metodă explicită în patru pași cu 4 ordine de precizie.

Diagrama bloc a procedurii de rezolvare a unei ecuații diferențiale prin metoda Runge-Kutta este prezentată în Figura 6.

F(x , y) - funcție dată - trebuie

fi descris separat.

Parametrii de intrare:
X0, X K - inițială și finală

independent

variabil;

Y 0 valoare y 0 din starea initiala

y (x 0 )= y 0 ;

N este numărul de segmente de partiție;

Parametri de ieșire:

Y - matrice de valori ale soluției dorite

la nodurile grilei.

2. 6. Rezolvarea problemei enunțate prin metodele Euler și Runge-Kutta de ordinul al 4-lea

2. 6. 1. Metoda Euler

1. Construim axele de coordonate;

2. Mark A (1,2; 1) primul punct al curbei integrale;

4. Construirea unei tangente l 0 în punctul A la un unghi α 0 ;

5. Găsiți x 1 cu formula: x i \u003d x 0 + ih, unde h etapa de integrare

x 1 \u003d 1,2 + 1 0,1 \u003d 1,3

6. Desenați o linie dreaptă x \u003d x 1 \u003d 0,1 până la intersecția cu linia l 0 , marcați punctul B (x 1 ; y 1 );

7. Căutând punctul y B :

Dintr-un triunghi dreptunghic ABC,

Δ y = y 1 y 0 ,

Δx \u003d x 1 x 0 \u003d h,

f (x 0 ; y 0 ) = (y 1 y 0 )/ h =>

y 1 \u003d y 0 + h (f (x 0; y 0 )) \u003d 1 + 0,1 f (1,2; 1) \u003d 1-0,545454 \u003d 0,945

De aici și ideea B are coordonatele (1,3; 0,945).

2. 6 .2. Metoda Runge-Kutta 4 ordine

1. Construim axele de coordonate;

2. Marcați A(1,2; 1) primul punct al curbei integrale;

3. Căutăm unghiul de înclinare al tangentei la grafic în punct A :

4. Construirea unei tangente l 0 în punctul A la un unghi α 0 ;

5. Găsiți x 1 după formula: x i \u003d x 0 + ih

x 1 \u003d 1, 2 + 1 0, 1 \u003d 1, 3;

  1. Găsim prin formulele:

k1=0,1 f(1,2;1)=0,1 (-0,545454)= - 0,0545

k2=0,1 f(1,2+0,1/2;1-0,0545/2)= -0,0540

K3=0,1 f(1,2+0,1/2;1-0,0540/2)= - 0,0540

K4=0,1 f(1,2+0,1;1-0,0540)= - 0,0535

∆y 1 =(- 0,0545+2 (-0,0540)+2 (-0,0540) - 0,0535)/6= - 0,054

∆y 2 \u003d 1- 0,054 \u003d 0,946

Prin urmare, următorul punct al graficului de decizie are coordonate (1,3; 0,946)

3. Algoritm pentru rezolvarea problemei

3.1.Algoritmi de subrutine

3.1.1 Rutina Euler

3.1.2 Subrutina metodei de ordine Runge-Kutta 4

3. 1. 3. Subprogramul soluţiei generale

3. 2. Algoritmul funcției


3. 3. Algoritmul programului

4. Forma programului

  1. Lista de programe

Dim j() As Single

Dim x() As Single

Dim y() Ca Single

Dim o() As Single

Privat n, i Ca întreg

Privat xk, x0, kx, ky Ca single

Privat k, k1, k2, k3, k4 Ca Single

Privat h, max, min, y0 Ca Single

Funcția privată f(a, b Ca Single) Ca Single

f = -a * b / (a ​​+ 1)

funcția finală

Funcția privată f1(x Ca Single) Ca Single

funcția finală

Sub Eiler privat()

ReDim x(n)

ReDim j(n)

j(0) = y0

Pentru i = 0 La n

x(i) = x0 + h * i

Apoi eu

Pentru i = 0 până la n - 1

Apoi eu

end sub

Curge secundar privat()

ReDim y(n)

y(0) = y0

Pentru i = 0 La n

x(i) = x0 + h * i

Apoi eu

Pentru i = 0 până la n - 1

k1 = h * f(x(i), y(i))

y(i + 1) = y(i) + k

MSFlexGrid1.TextMatrix(1, 2) = Str(y0)

Apoi eu

end sub

Sub Obhee privat()

ReDim o(n)

Pentru i = 0 La n

o(0) = y0

x(i) = x0 + h * i

o(i) = f1(x(i))

Apoi eu

end sub

Subcomandă privată1_Click()

x0 = Val(Text1.Text)

xk = Val(Text2.Text)

h = Val(Text4.Text)

y0 = Val(Text3.Text)

n = (xk - x0) / h

Label6.Caption = Str(x0)

Label5.Caption = Str(xk)

MSFlexGrid1.Rânduri = n + 2

MSFlexGrid1.TextMatrix(0, 1) = „Eleler”

MSFlexGrid1.TextMatrix(0, 2) = „Runge-Kutt”

MSFlexGrid1.TextMatrix(0, 3) = „Soluție generală”

Eiler

Runge

Obhee

max = y0

min = y0

Pentru i = 0 La n

Dacă j(i) > max Atunci

max = j(i)

Încheiați dacă

Dacă j(i)< min Then

min = j(i)

Încheiați dacă

Dacă y(i) > max Atunci

max = y(i)

Încheiați dacă

Dacă y(i)< min Then

min = y(i)

Încheiați dacă

Dacă o(i) > max Atunci

max = o(i)

Încheiați dacă

Dacă o(i)< min Then

min = o(i)

Încheiați dacă

Apoi eu

Label4.Caption = Str(max)

Label7.Caption = Str(min)

Poza1.Cls

Pentru i = 1 La n - 1

X1 = 720 + Rotunzi(kx * (x(i - 1) - x0))

X2 = 720 + rotund(kx * (x(i) - x0))

X1 = 720 + Rotunzi(kx * (x(i - 1) - x0))

X2 = 720 + rotund(kx * (x(i) - x0))

Apoi eu

end sub

Subcomandă privată2_Click()

end sub

  1. Rezolvarea problemelor în MathCad

Euler

Runge Kutt

General

Concluzie

Dintre cele două metode (Euler și Runge-Kutta), conform rezultatelor obținute, metoda Runge-Kutta s-a dovedit a fi mai precisă (comparând cu soluția generală). Acest lucru se explică prin faptul că, spre deosebire de metoda Euler din metoda Runge-Kutta, pasul nu este împărțit în 4 segmente, în urma cărora eroarea metodei devine mai mică.

La finalizarea cursului, am finalizat toate sarcinile: am scris un program pentru rezolvarea acestei ecuații diferențiale prin două metode numerice în program Visual Basic , a testat soluția cu aplicația MathCad , a comparat rezultatele obținute prin diferite metode cu soluția generală. Cred că am atins pe deplin scopul acestui curs.


tg(α) = f(x,y)

Y i+1 = Y i + h ∙ F(x, Y i )

x = X0 + i ∙ h

i = 0, …, N - 1

h = (Xk X0)/N

Sfârşit

MSFlexGrid1.TextMatrix(i + 2, 2) = Str(y(i + 1))

MSFlexGrid1.TextMatrix(1, 2) = Str( y(0))

kx = (6600 - 720) / (xk - x0)

ky = (6600 - 1120) / (max - min)

Label4.Caption = Str(max)

Label7.Caption = Str(min)

o(i)=min

o(i)

MSFlexGrid1.TextMatrix(i + 1, 3) = Str(o(i))

y(i)=max

Eiler(X0, Xk, Y0, N, Y)

x i +1

x i

yi+1

y=y(x)

k2=h*F(x+h/2, Y i +k1/2)

K1=h*F(x,Y i )

y(i) > max

j(i)=min

x = X0 + i ∙ h

i = 0, …, N-1

h = (Xk X0)/N

Rynge4(X0, Xk, Y0, N, Y)

Comanda 2

o(i)=max

o(i) > max

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

k4= h*F(x+h, Yi +k3)

k3= h*F(x+h/2, Yi +k2/2)

Yi+1 = Yi+k

k1 = h * f(x(i), y(i))

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

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

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

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

y(i + 1) = y(i) + k

i = 0, …, n-1

x(i) = x0 + h * i

i = 0, …, n

ReDim y(n)

g(0) = y0

start

y0,x0,xk,h

n = (xk - x0) / h

max = y0

min = y0

j(i) > max

i = 0, … n

Eiler

Runge

Obshee

MSFlexGrid1.Rânduri = n + 2

MSFlexGrid1.TextMatrix(0, 0) = „x”

MSFlexGrid1.TextMatrix(0, 1) = "Euler "

MSFlexGrid1.TextMatrix(0, 2) = "Runge- Kutta"

MSFlexGrid1.TextMatrix(0, 3) = "Decizie comună"

Label6.Caption = Str(x0)

Label5.Caption = Str(xk)

j(i)=max

X1 = 720 + Rotunzi(kx * (x(i - 1) - x0))

X2 = 720 + rotund(kx * (x(i) - x0))

Y1 = 6600 - Rotund(ky * (o(i - 1) - min))

Y2 = 6600 - Rotund(ky * (o(i) - min))

Sfârşit

j(i)

X1 = 720 + Rotunzi(kx * (x(i - 1) - x0))

X2 = 720 + rotund(kx * (x(i) - x0))

Y1 = 6600 - Rotund(ky * (y(i - 1) - min))

Y2 = 6600 - rotund(ky * (y(i) - min))

X1 = 720 + Rotunzi(kx * (x(i - 1) - x0))

X2 = 720 + rotund(kx * (x(i) - x0))

Y1 = 6600 - Rotund(ky * (j(i - 1) - min))

Y2 = 6600 - Rotund(ky * (j(i) - min))

i =1 , …, n-1

Sfârşit

f1 = y0 / ((x0 + 1) * Exp(-x0)) * (x + 1) * Exp(-x)

f1(x)

MSFlexGrid1.TextMatrix(1, 0) = Str(x0)

MSFlexGrid1.TextMatrix(i + 2, 0) = Str(x(i + 1))

MSFlexGrid1.TextMatrix(i + 2, 1) = Str(j(i + 1))

MSFlexGrid1.TextMatrix(1, 1) = Str(j(0))

j(i + 1) = j(i) + h * f(x(i), j(i))

i = 0, …, n-1

x(i) = x0 + h * i

i = 0, …, n

ReDim x(n)

ReDim j(n)

j(0) = y0

Eiler

Sfârşit

x(i) = x0 + h * i

o(i) =f1(x(i))

i = 0, …, n

ReDimo(n)

o(0) = y0

Obhee

Sfârşit

MSFlexGrid16

Imaginea 1

f = -a * b / (a ​​+ 1)

f(a,b)

Runge

Labe71

Text2

Text1

Labe41

Labe31

Eticheta 1

Text3

Labe21

Eticheta 6

Text4

Comanda 1

Eticheta4

Labe51

Eticheta 91

Eticheta 10

Sfârşit

Imagine1.Linie (X1, Y1)-(X2, Y2), RGB(0, 200, 0)

Imagine 1.Linie (X1, Y1)-(X2, Y2), RGB(500, 70, 90)

Imagine 1.Linie (X1, Y1)-(X2, Y2), RGB(400, 100, 12)

y(i)

y(i)=min

Să fie dată o ecuație diferențială de ordinul întâi

cu starea initiala

y(x0) = y0.

Alegem un pas h și introducem notația:

x i \u003d x 0 + i. h și y i = y(x i),

unde i = 0, 1, 2, …

x i – noduri de grilă,

y i - valoarea funcţiei integrale în noduri .

Similar metodelor descrise mai sus, ecuația diferențială este rezolvată. Diferența constă în împărțirea pasului în 4 părți.

Conform metodei Runge-Kutta de ordinul al patrulea, valori succesive y i funcția dorită y sunt determinate de formula:

, i = 0, 1, 2, …

si numerele k 1 (i), k 2 (i), k 3 (i), k 4 (i) la fiecare pas sunt calculate prin formulele:

Aceasta este o metodă explicită în patru etape de ordinul al patrulea de precizie.

Metodele Runge-Kutta sunt ușor de programat și au o precizie și stabilitate semnificative pentru o gamă largă de probleme.

Figura 6 este o diagramă a procedurii RUNGE(X0, XK, Y0, N, Y) pentru a rezolva problema Cauchy prin metoda Runge-Kutta descrisă mai sus.


i = 0, …, N-1

K1 = h * F(x, Yi)

K2 = h * F(x + h/2, Yi + K1 / 2)

K3 = h * F(x + h/2, Yi + K2 / 2)

K4 = h * F(x + h, Yi + K3)

K = (K1 + 2*K2 + 2*K3 + K4) / 6

Figura 6 - Organigrama procedurii RUNGE

Figura 7 prezintă o diagramă bloc a algoritmului programului principal pentru rezolvarea problemei Cauchy și obținerea de rezultate cu un număr fix de segmente de partiție N. În programul principal, procedura se numește RUNGE(X0, XK, Y0, N, Y), calcularea valorii funcției dorite y j la puncte xj prin metoda Runge-Kutta.

Datele inițialeîn această sarcină sunt:

X0, XK – valorile inițiale și finale ale variabilei independente;

Y0 - valoare y 0 din starea initiala y(x0) = y0;

N este numărul de segmente de partiție.

Rezultatele programului afisat in doua coloane:

X este o matrice de valori ale nodurilor grilei;

Y este o matrice de valori ale soluției dorite la nodurile grilei corespunzătoare.


Introduceți X0, XK, Y0, N


RUNGE(X0,XK,Y0,N,Y)


i=0…N


Ieșire X, Y i


Figura 7 - Diagrama bloc a algoritmului programului principal pentru rezolvarea problemei Cauchy cu un număr fix de segmente de partiție N


Rezolvarea ecuațiilor diferențiale în MathCad

Figura 8 - Un exemplu de rezolvare a unei ecuații diferențiale prin metoda

Runge-Kutta 4 comenzi în mediul MathCad.

Funcții de trasare în mediul Visual Basic

Sarcina de a tabula funcții și de a construi grafice ale acestora este una dintre sarcinile principale în procesul de rezolvare a ecuațiilor diferențiale obișnuite. Să luăm în considerare această problemă mai detaliat.

Sarcină.

Trasează funcția y=sin(x) pe segment. Faceți pasul de tabulare egal cu h.

Soluţie.

Pentru a reprezenta graficul unei funcții în mediul Visual Basic, este convenabil să folosiți unele componente grafice.

Figura 9 - Locația componentelor principale în fereastra General



Componenta Picture Box() este folosită ca container pentru trasare. Este o matrice de puncte (pixeli) și este posibil să controlați culoarea fiecărui punct individual. Coordonatele oricărui punct sunt determinate de o pereche de numere întregi - numărul său de serie în linia X și numărul de serie al liniei din interiorul obiectului Y. Astfel, coordonatele colțului din stânga sus al componentei sunt (0, 0) . Numărul de puncte pe linie și numărul de linii sunt determinate de dimensiunea componentei.

Figura 10 - Coordonatele obiectului PictureBox

Pe fig. 10 arată locația axelor și coordonatele punctelor de colț ale obiectului.

Componenta Line() este utilizată pentru a reprezenta axele și segmentele graficului poliliniei al funcției.

Esența trasării este că funcția trebuie prezentată sub forma unui tabel (tabulat), apoi marcați punctele pe șablonul de diagramă și conectați-le între ele.

Algoritmul pentru construirea unui grafic de funcții este prezentat în Figura 12. Algoritmul poate fi modificat. În special, unele proceduri pot fi fuzionate, iar procedura poate fi schimbată în unele cazuri.

Să luăm în considerare algoritmul mai detaliat.

Înainte de implementarea algoritmului, este necesar să se descrie o funcție-subrutină pentru trasarea unui grafic. Acest lucru este necesar pentru a facilita modificarea programului. Dacă trebuie să reprezentați o altă funcție, va fi suficient doar să schimbați subrutina.

De asemenea, înainte de a trasa un grafic, trebuie să creați și să editați un formular. Un exemplu de dezvoltare a unui formular este prezentat în Figura 11. Pe formular, trebuie să plasați componente pentru introducerea datelor inițiale, o componentă pentru tipărirea unui tabel, un buton de comandă și un container pentru plasarea unei diagrame (PictureBox). În interiorul PictureBox, trebuie să desenați axele de coordonate folosind linii drepte și să plasați etichete pentru a înregistra limitele segmentului de valori ale argumentului funcției și extremele funcției pe segmentul în cauză.

Introducerea datelor inițiale se realizează în programul în cauză prin apăsarea butonului de comandă. Foarte des, introducerea datelor este implementată folosind componenta TextBox.

Este convenabil să se efectueze procedura de tabelare a funcției într-o buclă cu un parametru, deoarece numărul de puncte de grafic care trebuie calculat este cunoscut. Înainte de a executa procedura, trebuie să specificați numărul de rânduri din tabel.

Numărul de linii este calculat prin formula k=n+2, unde k este numărul de linii și n este numărul de tabulaturi. Numărul de linii trebuie să fie mai mare decât numărul de segmente cu 2, deoarece este necesar să se țină cont de punctul de pornire (zero) și de linia pentru scrierea titlurilor coloanei paginii.

În procedura de tabulare în sine, pot fi combinate două acțiuni - tabularea și calcularea extremelor. Această soluție este prezentată în lista de programe din Figura 13.

Principala dificultate în trasare este trecerea de la valoarea matematică a funcției și a argumentului la coordonatele ecranului utilizate pentru a reprezenta graficul. La rezolvarea acestei probleme este necesar să se țină cont de direcția opusă a axelor de pe graficul matematic și de pe obiectul PictureBox și de necesitatea de a scala imaginea.

Factorii de scalare a graficului sunt calculați folosind următoarele formule:

unde kx este factorul de scalare de-a lungul axei OX,

NPX este numărul de pixeli ai obiectului PictureBox alocați pentru trasarea graficului pe orizontală,

a este valoarea inițială a argumentului de segment al funcției,

b este valoarea finală a segmentului argumentului funcției.

,

unde Ky este factorul de scalare de-a lungul axei OY,

NPY - numărul de pixeli ai obiectului PictureBox alocați pentru trasarea verticală,

min este valoarea minimă a funcției,

max este valoarea maximă a funcției.

Translația coordonatelor matematice ale punctului curent în coordonatele ecranului se realizează conform formulelor:

zx = rotund(ox + (x(i) - a) * kx),

zy = Rotund(oy - (y(i) - Min) * ky),

unde zx, zy sunt coordonatele ecranului punctului curent,

ox, oy - coordonatele punctului de intersecție al axelor din componenta pictureBox,

x(i), y(i) – coordonatele matematice ale punctului curent,

kx, ky sunt factori de scalare.

În formula de calcul a coordonatei ecranului ordonatei punctului curent, semnul minus este utilizat pentru a lua în considerare direcția opusă a axelor (pe ecran și pe grafic).

Lista programului pentru construirea unui grafic al unei funcții este prezentată în Figura 13.

Exemple de formulare cu rezultatele programului pentru diferite date inițiale sunt prezentate în figurile 14 și 15.

Figura 11 - Un exemplu de dezvoltare a unui formular

Figura 12 - Algoritm pentru construirea unui grafic al unei funcții

Rem Descrierea variabilelor

Dim x() ca single, y() ca single

Privat ca singur

  • Serghei Savenkov

    un fel de recenzie „rare”... parcă s-ar grăbi undeva