\(\displaystyle{ \begin{cases} x'(t) = -3x + 4y, \ \ x(0)=1, \\ y'(t) = -2x +3y, \ \ y(0) = 2, \ \ 0\leq t < 1 \end{cases} }\)
Kod: Zaznacz cały
function sys_rk4(f,g,a,b,x0,y0,n)
% Rozwiązanie układ równań różniczkowych x'=f(t,x,y), x(a)=x0 i
% y'=g(t,x,y), y(a)=y0 metodą Runge-Kutta rzędu 4.
fprintf('\n')
disp(' Rozwiązanie układu równań różniczkowych metodą Runge-Kutta rzędu4')
h=(b-a)/n;
x=x0;
y=y0;
disp('_____________________________________________________________________')
disp(' t x y x(t) y(t) |x-x(t)| |y-y(t)| ')
disp('_____________________________________________________________________')
fprintf('\n')
fprintf('%8.2f %10.3f %10.3f %10.3f %10.3f %8.2f %8.2f\n', a, x,y, x, y, 0, 0)
for i=1:n
t=a+(i-1)*h;
k1=feval(f,t,x,y);
c1=feval(g,t,x,y);
k2=feval(f,t+h/2,x+h*k1/2,y+h*c1/2);
c2=feval(g,t+h/2,x+h*k1/2,y+h*c1/2);
k3=feval(f,t+h/2,x+h*k2/2,y+h*c2/2);
c3=feval(g,t+h/2,x+h*k2/2,y+h*c2/2);
k4=feval(f,t+h,x+h*k3,y+h*c3);
c4=feval(g,t+h,x+h*k3,y+h*c3);
x=x+h*(k1+2*k2+2*k3+k4)/6;
y=y+h*(c1+2*c2+2*c3+c4)/6;
t=t+h;
xs='n';
ys='n';
if (xs=='n' & ys=='n')
fprintf('%8.2f %10.3f %10.3f\n',t,x,y)
elseif (ys~='n' & xs=='n')
err2=abs(ys-y);
fprintf('%8.2f %10.3f %10.3f %10.3f %10.1e\n', t, x, y, ys, err2)
else
err1=abs(xs-x);
err2=abs(ys-y);
fprintf('%8.2f %10.3f %10.3f %10.3f %10.3f %10.1e %10.1e\n', t, x, y, xs, ys, err1, err2)
end
end
function f = f1(t,x,y)
f= -3^x +4^y;
function g = g1(t,x,y)
g = -2*x + 3*y;
Kod: Zaznacz cały
Rozwiązanie układu równań różniczkowych metodą Runge-Kutta rzędu 4
________________________________________________________________________________________________________________________________
t x y x(t) y(t) |x-x(t)| |y-y(t)|
________________________________________________________________________________________________________________________________
0.00 -1.0000 1.0000 -1.0000 1.0000 0.000 0.000
0.10 -0.7627 0.9094 0.9094 1.29e-005
0.20 -0.1418 0.8601 0.8602 3.61e-005
0.30 1.0815 0.9009 0.9010 7.50e-005
0.40 3.2231 1.1068 1.1070 1.38e-004
0.50 6.7345 1.5909 1.5911 2.35e-004
0.60 12.2569 2.5203 2.5207 3.85e-004
0.70 20.6936 4.1388 4.1394 6.09e-004
0.80 33.3114 6.7976 6.7985 9.42e-004
0.90 51.8779 10.9982 10.9986 1.43e-003
1.00 78.8502 17.4517 17.4539 2.13e-003
>> # Octave 7.1.0, Mon Aug 05 17:59:08 2024 GMT
Kod: Zaznacz cały
function sys_rk4(f,g,a,b,x0,y0,n)
% Rozwiązanie układ równań różniczkowych x'=f(t,x,y), x(a)=x0 i
% y'=g(t,x,y), y(a)=y0 metodą Runge-Kutta rzędu 4.
fprintf('\n')
disp(' Rozwiązanie układu równań różniczkowych metodą Runge-Kutta rzędu4')
h=(b-a)/n;
x=x0;
y=y0;
disp('_____________________________________________________________________')
disp(' t x y x(t) y(t) |x-x(t)| |y-y(t)| ')
disp('_____________________________________________________________________')
fprintf('\n')
fprintf('%8.2f %10.3f %10.3f %10.3f %10.3f %8.2f %8.2f\n', a, x,y, x, y, 0, 0)
for i=1:n
t=a+(i-1)*h;
k1=feval(f,t,x,y);
c1=feval(g,t,x,y);
k2=feval(f,t+h/2,x+h*k1/2,y+h*c1/2);
c2=feval(g,t+h/2,x+h*k1/2,y+h*c1/2);
k3=feval(f,t+h/2,x+h*k2/2,y+h*c2/2);
c3=feval(g,t+h/2,x+h*k2/2,y+h*c2/2);
k4=feval(f,t+h,x+h*k3,y+h*c3);
c4=feval(g,t+h,x+h*k3,y+h*c3);
x=x+h*(k1+2*k2+2*k3+k4)/6;
y=y+h*(c1+2*c2+2*c3+c4)/6;
t=t+h;
xs='n';
ys='n';
if (xs=='n' & ys=='n')
fprintf('%8.2f %10.3f %10.3f\n',t,x,y)
elseif (ys~='n' & xs=='n')
err2=abs(ys-y);
fprintf('%8.2f %10.3f %10.3f %10.3f %10.1e\n', t, x, y, ys, err2)
else
err1=abs(xs-x);
err2=abs(ys-y);
fprintf('%8.2f %10.3f %10.3f %10.3f %10.3f %10.1e %10.1e\n', t, x, y, xs, ys, err1, err2)
end
end
function f = f1(t,x,y)
f= -3^x +4^y;
function g = g1(t,x,y)
g = -2*x + 3*y;
Kod: Zaznacz cały
Rozwiązanie układu równań różniczkowych metodą Runge-Kutta rzędu 4
________________________________________________________________________________________________________________________________
t x y x(t) y(t) |x-x(t)| |y-y(t)|
________________________________________________________________________________________________________________________________
0.00 -1.0000 1.0000 -1.0000 1.0000 0.000 0.000
0.10 -0.7627 0.9094 0.9094 1.29e-005
0.20 -0.1418 0.8601 0.8602 3.61e-005
0.30 1.0815 0.9009 0.9010 7.50e-005
0.40 3.2231 1.1068 1.1070 1.38e-004
0.50 6.7345 1.5909 1.5911 2.35e-004
0.60 12.2569 2.5203 2.5207 3.85e-004
0.70 20.6936 4.1388 4.1394 6.09e-004
0.80 33.3114 6.7976 6.7985 9.42e-004
0.90 51.8779 10.9982 10.9986 1.43e-003
1.00 78.8502 17.4517 17.4539 2.13e-003
>> # Octave 7.1.0, Mon Aug 05 17:59:08 2024 GMT