ASSIGNMENT-2
RUNGE-KUTTA METHOD TO SOLVE ORDINARY DIFFERENTIAL EQUATIONS
PROBLEM STATEMENT:
y’ = -2xy2 ; y(0) = 1
Find y at x = 1 using Runge – Kutta method of order 4 by solving the given IVP. Take step
size of 0.2
MODEL:
Runge – Kutta method of order 4 gives value of the functions as follows:
k1 = h f (xn , yn)
k2 = h f(xn + h/2, yn + k1/2)
k3 = h f(xn + h/2, yn + k2/2)
k4 = h f(xn + h, yn + k3)
yn+1 = yn + (k1 + 2k2 +2k3 + k4)/6
The solution to this ODE after solving was y = 1/(x2 + 1) + C
Using the given initial conditions, we found that C = 0;
Thus, the equation became y = 1/(x2 + 1), and y at x = 1 = 0.5, which is the same solution we
found using the RK method.
CODE:
y_dash = @(x,y) -2*x*y^2;
x_initial = 0;
y_initial = 1;
step_length = 0.2;
i = 1;
while x_initial<=1
k1 = step_length*y_dash(x_initial,y_initial);
k2 = step_length*y_dash((x_initial+step_length/2),(y_initial + k1/2));
k3 = step_length*y_dash((x_initial+step_length/2),(y_initial + k2/2));
k4 = step_length*y_dash((x_initial+step_length),(y_initial + k3));
km = (k1+2*k2+2*k3+k4)/6;
x(i) = x_initial;
y(i) = y_initial;
y_initial = y_initial + km*step_length;
x_initial = x_initial + step_length;
i = i+1;
end
disp(y_initial)
RESULTS AND DISCUSSION:
The value of y at x=1 was found to be 0.8095
PROBLEM STATEMENT:
[y1’, y2’] = [ 9y1 + 24y2 + (5 cos(t) – (1/3) sin(t), -24y1 -51y2 -9 cos(t) + (1/3) sin(t)]
y0 = [4/3, 2/3]
Solve the system of differential equations by 4th Order RK method from t = 0 to t = 1 using
step size (h) of 0.05, and find the area under the curve
MODEL:
Runge – Kutta method of order 4 gives value of the functions as follows:
k1 = h f (xn , yn)
k2 = h f(xn + h/2, yn + k1/2)
k3 = h f(xn + h/2, yn + k2/2)
k4 = h f(xn + h, yn + k3)
yn+1 = yn + (k1 + 2k2 +2k3 + k4)/6
Here, since there are 3 variables in each equation (y1, y2, t), each call becomes
k1 = h f(yn1 , yn2, t)
l1 = h g(xn1 , yn2, t)
k2 = h f(yn1 + k1/2, yn2 + l1/2, t + h/2)
k2 = h g(yn1 + k1/2, yn2 + l1/2, t + h/2)
and so on.
Where f = y1’ , and g = y2’
CODE:
% Defining the system of differential equations
f1 = @(t, y1, y2) 9*y1 + 24*y2 + 5*cos(t) - (1/3)*sin(t);
f2 = @(t, y1, y2) -24*y1 - 51*y2 - 9*cos(t) + (1/3)*sin(t);
% Initialising the conditions
y1_0 = 4/3;
y2_0 = 2/3;
t0 = 0;
tf = 1;
h = 0.05;
% Number of steps
N = (tf - t0) / h;
% Initialize arrays to store the solutions
t = t0:h:tf;
y1 = zeros(1, N+1);
y2 = zeros(1, N+1);
% Set initial conditions
y1(1) = y1_0;
y2(1) = y2_0;
% RK method of fourth order
for i = 1:N
k1_y1 = h * f1(t(i), y1(i), y2(i));
k1_y2 = h * f2(t(i), y1(i), y2(i));
k2_y1 = h * f1(t(i) + h/2, y1(i) + k1_y1/2, y2(i) + k1_y2/2);
k2_y2 = h * f2(t(i) + h/2, y1(i) + k1_y1/2, y2(i) + k1_y2/2);
k3_y1 = h * f1(t(i) + h/2, y1(i) + k2_y1/2, y2(i) + k2_y2/2);
k3_y2 = h * f2(t(i) + h/2, y1(i) + k2_y1/2, y2(i) + k2_y2/2);
k4_y1 = h * f1(t(i) + h, y1(i) + k3_y1, y2(i) + k3_y2);
k4_y2 = h * f2(t(i) + h, y1(i) + k3_y1, y2(i) + k3_y2);
y1(i+1) = y1(i) + (k1_y1 + 2*k2_y1 + 2*k3_y1 + k4_y1) / 6;
y2(i+1) = y2(i) + (k1_y2 + 2*k2_y2 + 2*k3_y2 + k4_y2) / 6;
end
% Display the results
fprintf('t\t\ty1\t\ty2\n');
for i = 1:N+1
fprintf('%.2f\t%.4f\t%.4f\n', t(i), y1(i), y2(i));
end
% Plot the results
figure;
subplot(2,1,1);
plot(t, y1, '-o');
xlabel('t');
ylabel('y1');
title('Solution for y1');
RESULTS AND DISCUSSION:
Time (seconds) Y2 Y1
0.00 1.3333 0.6667
0.05 1.7364 -0.5578
0.10 1.7122 -0.8703
0.15 1.5727 -0.9029
0.20 1.4141 -0.8550
0.25 1.2644 -0.7888
0.30 1.1305 -0.7229
0.35 1.0126 -0.6623
0.40 0.9093 -0.6079
0.45 0.8186 -0.5593
0.50 0.7388 -0.5156
0.55 0.6682 -0.4762
0.60 0.6057 -0.4404
0.65 0.5499 -0.4076
0.70 0.4998 -0.3774
0.75 0.4547 -0.3492
0.80 0.4137 -0.3229
0.85 0.3761 -0.2980
0.90 0.3416 -0.2744
0.95 0.3096 -0.2517
1.00 0.2797 -0.2299
