Diferencijalne jednačine

Razmatraćemo Košijev problem

pri čemu ćemo određivati približne vrednosti tačnog rešenja u istaknutim tačkama

Diferecijalne.png

Na osnovu Njutn-Lajbnicove teoreme znamo da je

odakle je

Ako iskoristimo neku kvadraturnu formulu za aproksimaciju integrala na desnoj strani, tada možemo računati vrednost na osnovu vrednosti . Svakako ćemo uzeti da je U zavisnosti od odabira kvadraturne formule ćemo imati različite metode za dobijanje numeričkog rešenja datog Košijevog problema.

Ojlerova metoda

Približne vrednosti u čvorovima se računaju po formuli:

za Ovde smo koristili jednostavnu formulu pravougaonika za aproksimaciju odgovarajućeg integrala. U nastavku je data jednostavna implementacije ove metode. Druga kolona u izlaznom ispisu je kolona vrednosti dobijenih Ojlerovom metodom, a treća je kolona vrednosti tačnog rešenja u čvorovima.

% function Euler(f,a,b,h,y0,y)

% fplot(y,[a,b],'-r');

% hold on

% x = a:h:b;

% n = length(x);

% ind = 0:n-1;

% yNumerical = zeros(1,n);

% yNumerical(1) = y0;

% for i = 1:n-1

% yNumerical(i+1) = yNumerical(i)+h*f(x(i),yNumerical(i));

% end

% scatter(x,yNumerical,'b');

% legend('Tacno resenje','Numericko resenje','Location','NorthWest');

% hold off

% Res(:,1) = x;

% Res(:,2) = yNumerical;

% Res(:,3) = y(x);

% disp(Res)

% end

format short

f = @(x,y)x+y;

a = 0;

b = 1;

h = 0.1;

y0 = 1;

y = @(x)2*exp(x)-x-1;

Euler(f,a,b,h,y0,y);

         0    1.0000    1.0000
    0.1000    1.1000    1.1103
    0.2000    1.2200    1.2428
    0.3000    1.3620    1.3997
    0.4000    1.5282    1.5836
    0.5000    1.7210    1.7974
    0.6000    1.9431    2.0442
    0.7000    2.1974    2.3275
    0.8000    2.4872    2.6511
    0.9000    2.8159    3.0192
    1.0000    3.1875    3.4366

Metoda Runge Kuta

Postoji više formula Runge-Kuta različite tačnosti. Ovde ćemo navesti dve, Runge-Kutovu formulu drugog i četvrtog reda.

Drugog reda

Ova metoda se može posmatrati kao vrlo jasna modifikacija Ojlerove metode:

Implementacija:

% function Runge2(f,a,b,h,y0,y)

% x = a:h:b;

% n = length(x);

% ind = 0:n-1;

% yNumerical = zeros(1,n);

% yNumerical(1) = y0;

% for i = 1:n-1

% yNumerical(i+1) = yNumerical(i)+h*f(x(i)+h/2,yNumerical(i)+(h/2)*f(x(i),yNumerical(i)));

% end

% Res(:,1) = x;

% Res(:,2) = yNumerical;

% Res(:,3) = y(x);

% disp(Res)

% end

format short

f = @(x,y)x+y;

a = 0;

b = 1;

h = 0.1;

y0 = 1;

y = @(x)2*exp(x)-x-1;

Runge2(f,a,b,h,y0,y);

         0    1.0000    1.0000
    0.1000    1.1100    1.1103
    0.2000    1.2421    1.2428
    0.3000    1.3985    1.3997
    0.4000    1.5818    1.5836
    0.5000    1.7949    1.7974
    0.6000    2.0409    2.0442
    0.7000    2.3231    2.3275
    0.8000    2.6456    2.6511
    0.9000    3.0124    3.0192
    1.0000    3.4282    3.4366

Četvrtog reda

Rešenja se u ovom slučaju dobijaju vezom

gde je redom

Za sve prethodne metode je moguće odrediti aproksimaciju greške preko formule

gde su i približne vrednosti kada računamo sa korakom i , redom, a parametar koji zavisi od konkretne metode. U slučaju Ojlerove metode, , a u slučaju metoda Runge-Kuta, , odnosno .

U nastavku je data implementacija metode Runge-Kuta reda 4 za koju je određena i aproksimacija greške. Primetimo da formula za ocenu greške ima smisla samo za svaki drugi čvor. Za ostale čvorove možemo tu aproksimaciju samo da procenimo na neki način, na primer, tako što uzmemo aritmetičku sredinu susednih grešaka. U zapisu će biti ispisane kolone rešnenja prvo sa korakom , a zatim sa korakom . Poslednje kolona će kao i ranije biti kolona sa vrednostima tačnog rešanja u čvorovima.

% function Runge4(f,a,b,h,y0,y)

% x = a:h:b;

% n = length(x);

% yNumerical = zeros(1,n);

% yNum2h = zeros(1,n);

% ApproxErr = zeros(1,n);

% yNumerical(1) = y0;

% yNum2h(1) = y0;

% i = 1;

% while i<=n-1

% K1 = h*f(x(i),yNumerical(i));

% K2 = h*f(x(i)+h/2,yNumerical(i)+K1/2);

% K3 = h*f(x(i)+h/2,yNumerical(i)+K2/2);

% K4 = h*f(x(i)+h,yNumerical(i)+K3);

% yNumerical(i+1) = yNumerical(i)+(1/6)*(K1+2*K2+2*K3+K4);

% if (mod(i,2)==0)

% k1 = (2*h)*f(x(i-1),yNum2h(i-1));

% k2 = (2*h)*f(x(i-1)+h,yNum2h(i-1)+k1/2);

% k3 = (2*h)*f(x(i-1)+h,yNum2h(i-1)+k2/2);

% k4 = (2*h)*f(x(i-1)+2*h,yNum2h(i-1)+k3);

% yNum2h(i+1) = yNum2h(i-1)+(1/6)*(k1+2*k2+2*k3+k4);

% ApproxErr(i+1) = (abs(yNumerical(i+1)-yNum2h(i+1)))/15;

% ApproxErr(i) = (ApproxErr(i-1)+ApproxErr(i+1))/2;

% end

% i = i+1;

% end

% Res(:,1) = x;

% Res(:,2) = yNumerical;

% Res(:,3) = yNum2h;

% Res(:,4) = y(x);

% Res(:,5) = ApproxErr;

% disp(Res)

% end

format long

f = @(x,y)x+y;

a = 0;

b = 1;

h = 0.1;

y0 = 1;

y = @(x)2*exp(x)-x-1;

Runge4(f,a,b,h,y0,y);

  Columns 1 through 3

                   0   1.000000000000000   1.000000000000000
   0.100000000000000   1.110341666666667                   0
   0.200000000000000   1.242805141701389   1.242800000000000
   0.300000000000000   1.399716994125076                   0
   0.400000000000000   1.583648480161372   1.583635920000000
   0.500000000000000   1.797441277193677                   0
   0.600000000000000   2.044235924183866   2.044212912688000
   0.700000000000000   2.327503253193554                   0
   0.800000000000000   2.651079126584631   2.651041651557123
   0.900000000000000   3.019202827560142                   0
   1.000000000000000   3.436559488270332   3.436502273211870

  Columns 4 through 5

   1.000000000000000                   0
   1.110341836151295   0.000000171390046
   1.242805516320340   0.000000342780093
   1.399717615152007   0.000000590062092
   1.583649395282541   0.000000837344091
   1.797442541400256   0.000001185721908
   2.044237600781018   0.000001534099724
   2.327505414940953   0.000002016217446
   2.651081856984936   0.000002498335167
   3.019206222313899   0.000003156336199
   3.436563656918091   0.000003814337231

MATLAB nudi i ugrađene funkcije za numeričko rešavanje diferencijalnih jednačina. Neke od njih su ode15(), ode23(), ode45()... Ovde ćemo razmotriti primer poziva funkcije ode45() koja dolazi iz familije Runge-Kutovih metoda.

format long

f = @(x,y)x+y;

X = 0:0.1:1;

y0 = 1;

[x,y] = ode45(f,X,y0);

disp([x,y]);

                   0   1.000000000000000
   0.100000000000000   1.110341836666667
   0.200000000000000   1.242805517459487
   0.300000000000000   1.399717617040434
   0.400000000000000   1.583649398065255
   0.500000000000000   1.797442545244475
   0.600000000000000   2.044237605879241
   0.700000000000000   2.327505421514429
   0.800000000000000   2.651081865287580
   0.900000000000000   3.019206232636722
   1.000000000000000   3.436563669594182