Runge Kutta 4th order

Hello,

Here is the task that i have to solve:

y1' = y2

y2'=f(x,y1,y2) with y1(0)=0 and y2(0)=y20

where f(x,y1,y2) = -axy2-y1, a=0.03 and y20 = 0.25

and here is my matlab code:

--

function main
h = pi/10;
x = 0 : h : 2*pi;
n = length(x);
y(1,1)=0;  y(2,1) = 0,25;
for k = 1: (n-1)
                k1 = h * rung(x(k), y(1:2,k)); 
                k2 = h * rung(x(k) + h/2, y(1:2,k) + k1/2);
                k3 = h * rung(x(k) + h/2, y(1:2,k) + k2/2);
                k4 = h * rung(x(k) + h, y(1:2,k) + k3);
                y(1:2,k+1)=y(1:2,k) + (k1 + 2*k2 + 2*k3 + k4)/6;
end
figure(1), plot(x, y(2,:)) %everything from second row
figure(2), plot(x, y(1,:), x, ??, '*')
function f=rung(x,y)
f(1,1) = y(2); 
f(2,1) = -0,03*x*y(2) – y(1);

So, my question is - Is this correct or not? And what would be the best option for h and x.

And how can i found and use numerical solution? Because i know that instead of the question marks here figure(2), plot(x, y(1,:), x, ??, '*') should be the exact numerical solution.

Answer 1

[Erivelton Gualter]

Your code seems to have some issues. Try this one:

h = pi/10;
x = 0 : h : 2*pi;
n = length(x);

y = zeros(2,n);
y(:,1) = [0, 0.25];
for k = 1: (n-1)
                k1 = h * rung(x(k), y(:,k)); 
                k2 = h * rung(x(k)+h/2, y(:,k) + k1/2);
                k3 = h * rung(x(k)+h/2, y(:,k) + k2/2);
                k4 = h * rung(x(k)+h, y(:,k) + k3);
                y(:,k+1) = y(:,k) + (k1 + 2*k2 + 2*k3 + k4)/6;
end

subplot(211); plot(x, y(1,:))
subplot(212); plot(x, y(2,:))

function f = rung(x,y)
    f(1,1) = y(2); 
    f(2,1) = -0.03*x*y(2) - y(1);
end

[Melda Harlova]

Thank you very much! It seems to be ok. But i did not understand why you usе this: y = zeros(2,n) and is this solution is numeric.

Thank a lot again. :)

Answer 2

[Meysam Mahooti]

https://www.mathworks.com/matlabcentral/fileexchange/55430-runge-kutta-4th-order?s_tid=srchtitle

https://www.mathworks.com/matlabcentral/fileexchange/61130-runge-kutta-fehlberg-rkf78?s_tid=srchtitle

Answer 3

[Abhinanda]

function [xRK, yRK] = runge_kutta(y0, xspan, h) xRK = xspan(1):h:xspan(2); yRK = zeros(size(xRK)); yRK(1) = y0;

    for i = 1:(length(xRK)-1)
        k1 = h * (-2*xRK(i)^3 + 12*xRK(i)^2 - 20*xRK(i) + 8.5);
        k2 = h * (-2*(xRK(i)+h/2)^3 + 12*(xRK(i)+h/2)^2 - 20*(xRK(i)+h/2) + 8.5);
        k3 = h * (-2*(xRK(i)+h/2)^3 + 12*(xRK(i)+h/2)^2 - 20*(xRK(i)+h/2) + 8.5);
        k4 = h * (-2*(xRK(i)+h)^3 + 12*(xRK(i)+h)^2 - 20*(xRK(i)+h) + 8.5);
        yRK(i+1) = yRK(i) + (k1 + 2*k2 + 2*k3 + k4) / 6;
    end
end

% Run Runge-Kutta method [xRK, yRK] = runge_kutta(10, [0, 4], 0.5); plot(xRK, yRK, '-o'); legend('Analytical Solution', 'Euler Explicit', 'Runge-Kutta');

Answer 4

[Abhinanda]

function [xI, yEuler2] = euler_implicit(y0, xspan, h) xI = xspan(1):h:xspan(2); yEuler2 = zeros(size(xI)); yEuler2(1) = y0;

    for i = 1:(length(xI)-1)
        % Solve for yEuler2(i+1) using fsolve
        y_guess = yEuler2(i);  % Initial guess for fsolve
        f = @(y) y - yEuler2(i) - h * (-2*xI(i+1)^3 + 12*xI(i+1)^2 - 20*xI(i+1) + 8.5);
        yEuler2(i+1) = fsolve(f, y_guess);
    end
end

% Run Euler's implicit method [xI, yEuler2] = euler_implicit(10, [0, 4], 0.5); plot(xI, yEuler2, '-o'); legend('Analytical Solution', 'Euler Explicit', 'Runge-Kutta', 'Euler Implicit');

Answer 5

[Abhinanda]

% Define the ODE as a function handle odefun = @(x, y) -2*x^3 + 12*x^2 - 20*x + 8.5;

% Use ode23 to solve the ODE [xMATLAB, yMATLAB] = ode23(odefun, [0 4], 10);

% Plot all solutions together for comparison plot(xMATLAB, yMATLAB, '-o'); legend('Analytical Solution', 'Euler Explicit', 'Runge-Kutta', 'Euler Implicit', 'MATLAB ode23'); hold off;

Answer 6

[Abhinanda]

function [xI, yEuler2] = euler_implicit_manual(y0, xspan, h) xI = xspan(1):h:xspan(2); yEuler2 = zeros(size(xI)); yEuler2(1) = y0;

    for i = 1:(length(xI)-1)
        % Use an iterative approach to estimate the next y value
        y_current = yEuler2(i);
        x_next = xI(i+1);
        y_next_guess = y_current;  % Initial guess for next y
        tol = 1e-6;  % Tolerance level for convergence

        % Iterate until convergence
        while true
            % Implicit Euler formula: y_new = y_old + h * f(x_new, y_new)
            y_new = y_current + h * (-2*x_next^3 + 12*x_next^2 - 20*x_next + 8.5);
            if abs(y_new - y_next_guess) < tol
                break;
            end
            y_next_guess = y_new;
        end
        yEuler2(i+1) = y_new;
    end
end

% Run Euler's implicit method using the manual iteration [xI, yEuler2] = euler_implicit_manual(10, [0, 4], 0.5); plot(xI, yEuler2, '-o'); legend('Analytical Solution', 'Euler Explicit', 'Runge-Kutta', 'Euler Implicit (Manual)');

Answer 7

[Abhinanda]

% Define the ODE as a function handle odefun = @(x, y) -2*x^3 + 12*x^2 - 20*x + 8.5;

% Use ode45 to solve the ODE [xMATLAB, yMATLAB] = ode45(odefun, [0 4], 10);

% Plot all solutions together for comparison plot(xMATLAB, yMATLAB, '-o'); legend('Analytical Solution', 'Euler Explicit', 'Runge-Kutta', 'Euler Implicit (Manual)', 'MATLAB ode45'); hold off;

Answer 8

[Abhinanda]

% Define the rate constants k1 = 2.1; % L/(mol.s) k2 = 1.5; % L/(mol.s) k3 = 1.3; % L/(mol.s)

% Initial concentrations x0 = 1.0; % mol/L y0 = 0.2; % mol/L z0 = 0.0; % mol/L

% Time span for the simulation tspan = [0 3]; % seconds

% Define the system of ODEs as a function handle function dydt = reaction_system(t, y) x = y(1); y = y(2); z = y(3);

    dxdt = -k1*x + k2*y*z;
    dydt = k1*x - k2*y*z - 2*k3*y^2;
    dzdt = k3*y^2;

    dydt = [dxdt; dydt; dzdt];
end

% Solve the system of ODEs using ode45 [t, sol] = ode45(@reaction_system, tspan, [x0, y0, z0]);

% Extract the solutions for x, y, and z xVal = sol(:,1); yVal = sol(:,2); zVal = sol(:,3);

% Plot the concentrations over time figure; plot(t, xVal, '-o', 'DisplayName', 'X'); hold on; plot(t, yVal, '-x', 'DisplayName', 'Y'); plot(t, zVal, '-s', 'DisplayName', 'Z'); xlabel('Time (s)'); ylabel('Concentration (mol/L)'); title('Concentration of Species X, Y, Z over Time'); legend('show'); grid on;

Answer 9

[Abhinanda]

k1 = 2.1; k2 = 1.5; k3 = 1.3;

x0 = 1.0; y0 = 0.2; z0 = 0.0;

tspan = [0 3];

function dydt = reaction_system(t, y) x = y(1); y = y(2); z = y(3);

    dxdt = -k1*x + k2*y*z;
    dydt = k1*x - k2*y*z - 2*k3*y^2;
    dzdt = k3*y^2;

    dydt = [dxdt; dydt; dzdt];
end

[t, sol] = ode45(@reaction_system, tspan, [x0, y0, z0]);

xVal = sol(:,1); yVal = sol(:,2); zVal = sol(:,3);

figure; plot(t, xVal, '-o', 'DisplayName', 'X'); hold on; plot(t, yVal, '-x', 'DisplayName', 'Y'); plot(t, zVal, '-s', 'DisplayName', 'Z'); xlabel('Time (s)'); ylabel('Concentration (mol/L)'); title('Concentration of Species X, Y, Z over Time'); legend('show'); grid on;

Answer 10

[Abhinanda]

k1 = 10^-2; k2 = 10^4; k3 = 10^2;

x0 = 1.0; y0 = 0.2; z0 = 0.0;

tspan = [0 1000];

function dydt = reaction_system(t, y) x = y(1); y = y(2); z = y(3);

    dxdt = -k1*x + k2*y*z;
    dydt = k1*x - k2*y*z - 2*k3*y^2;
    dzdt = k3*y^2;

    dydt = [dxdt; dydt; dzdt];
end

[t1, sol1] = ode45(@reaction_system, tspan, [x0, y0, z0]); [t2, sol2] = ode15s(@reaction_system, tspan, [x0, y0, z0]);

xVal1 = sol1(:,1); yVal1 = sol1(:,2); zVal1 = sol1(:,3);

xVal2 = sol2(:,1); yVal2 = sol2(:,2); zVal2 = sol2(:,3);

figure; subplot(3,1,1); plot(t1, xVal1, '-o', 'DisplayName', 'X - ode45'); hold on; plot(t2, xVal2, '-x', 'DisplayName', 'X - ode15s'); xlabel('Time (s)'); ylabel('X Concentration (mol/L)'); legend('show');

subplot(3,1,2); plot(t1, yVal1, '-o', 'DisplayName', 'Y - ode45'); hold on; plot(t2, yVal2, '-x', 'DisplayName', 'Y - ode15s'); xlabel('Time (s)'); ylabel('Y Concentration (mol/L)'); legend('show');

subplot(3,1,3); plot(t1, zVal1, '-o', 'DisplayName', 'Z - ode45'); hold on; plot(t2, zVal2, '-x', 'DisplayName', 'Z - ode15s'); xlabel('Time (s)'); ylabel('Z Concentration (mol/L)'); legend('show');

Answer 11

[Abhinanda]

function dydx = second_order_ode(x, y) dydx = [y(2); 5 * y(1) - 6 * y(2)]; end

x_span = 0:0.1:2; y0 = [1; 2];

[x, ySol] = ode45(@second_order_ode, x_span, y0);

yVal1 = ySol(:, 1); y1Val = ySol(:, 2);

yAct = exp(2 * x);

figure; plot(x, yVal1, 'o-', 'DisplayName', 'Numerical Solution (y)'); hold on; plot(x, yAct, 'x-', 'DisplayName', 'Analytical Solution (y=exp(2x))'); xlabel('x'); ylabel('y'); legend('show'); grid on;

Answer 12

[Abhinanda]

function dydx = ode_system(x, y) dydx = [y(2); 5*y(2) - 6*y(1)]; end

x_vals = [0, 0.2, 0.3, 0.45, 0.57, 0.7, 0.81, 0.9, 0.96, 1]; y0 = [1; 2]; [xSol, ySol] = ode45(@ode_system, x_vals, y0);

yVal = ySol(:, 1); y1Val = ySol(:, 2); yAct = exp(2*x_vals);

plot(x_vals, yVal, '-o', x_vals, yAct, '--x'); legend('yVal', 'yAct');

Answer 13

[Abhinanda]

1 x spn [0.0.2,0.3,0.45,0.57,0.7.0.81,0.9,8.96,1];

2 [x,y] ode45(@vdp,x_span, [12])

yvaly(1:10,1); 4ylValy(1:10,2);

5 yAct 11: for 11:10

temp exp(2*x_span(1,1)); yAct [yAct, temp];

end

10 yAct yAct

11 plot(x,yval, "Color", "r")

12 hold on

13 plot(x, yAct, LineWidth-1,color="b")

14 hold off

15

16 function dydt vdp(x,y)

17 dydt [y(2); 5*y(2)-6*y(1)];

18 end