Runge Kutta 4th order

1 次查看（过去 30 天）
Melda Harlova
Melda Harlova 2019-5-8
回答： Abhinanda 2024-10-23
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.

请先登录，再回答此问题。

采纳的回答

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
  1 个评论
Melda Harlova
Melda Harlova 2019-5-9
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. :)

请先登录，再进行评论。

更多回答（11 个）

Abhinanda
Abhinanda 2024-10-23

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');

Abhinanda
Abhinanda 2024-10-23

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');

Abhinanda
Abhinanda 2024-10-23
% 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;
Abhinanda
Abhinanda 2024-10-23

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)');

Abhinanda
Abhinanda 2024-10-23
% 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;
Abhinanda
Abhinanda 2024-10-23

% 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;

Abhinanda
Abhinanda 2024-10-23

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;

Abhinanda
Abhinanda 2024-10-23

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');

Abhinanda
Abhinanda 2024-10-23
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;
Abhinanda
Abhinanda 2024-10-23
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');
Abhinanda
Abhinanda 2024-10-23
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

类别

在 Help Center 和 File Exchange 中查找有关 Introduction to Installation and Licensing 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by