Vectorizing a Large Function with a Command
Show older comments
Hello, everyone. I am trying to run this program, however, I have run into a big problem. As you can see below, the code is used to solve a 2DOF system of differential equations for the dissociation reaction of CO2. However, the differential equations are much too large to tediously sit through and change each operation in order to allow for solving the system as an array. I am wondering, there has to be a way to automate this process. To automatically "Vectorize" a function so it can be solved with matrices. In engineering, there are countless complex equations that would take considerable effort to change each operation to suite the array solver needs. I have also been running into a strange issue where the command window gives a "Too many input arguments" error for the function of dy1, when the variables chosen for the runge kutta integration for loop are all in the function itself.
I am very very close to solving this reaction, and doing so would help me move further along in my internal combustion engine research. I am open to any suggestions, and any help would be greatly appreciated. Thanks.
P
Code:
% Chemical Equilibrium
% Solve the following chemical equation:
% aCO2 -> n1[CO2]+n2[O2]+n3[CO]+n4[O]
% Goal: We want to be able to solve the chemical reaction above for all
% products assuming chemical equilibrium
% Assume a fraction of CO2 and O2 dissociates. We will call these parameters
% "y1" and "y2"
% Write out the mass balance:
% a -> a(1-y1)[CO2]+(1/2)ay1(1-y2)[O2]+ay1[CO]+ay1[CO]+ay1y2[O]
% From this, we can derive the total number of products:
% np = a[1-y1+(1/2)y1(1-y2)+y1+y1y2]
% np = (1/2)a[2+y1(1+y2)]
% From this, we define each mole fraction
% x1 = 2(1-y1)/(2+y1(1+y2))
% x2 = y1(1-y2)/(2+y1(1+y2))
% x3 = 2y1/(2+y1(1+y2))
% x4 = 2y1y2/(2+y1(1+y2))
% From here, we can plug these values into the equilibrium constant
% equation and solve for the time derivate of y. To solve for the
% differential equation, initialization parameters xij are used.
clearvars
syms K1 dK1 K2 dK2 y1 dy1 y2 dy2 P dP
% Initialization Parameters: These mole fraction and dissociation constant
% equations are being used to solve for the time derivatives of y1 and y2
x1j = 2*(1-y1)/(2+y1*(1+y2));
x2j = y1*(1-y2)/(2+y1*(1+y2));
x3j = y1*(1-y2)/(2+y1*(1+y2));
x4j = y1*(1-y2)/(2+y1*(1+y2));
% Differential Initialization Parameters
dx1j = diff(x1j,y1)*dy1+diff(x1j,y2)*dy2;
dx2j = diff(x2j,y1)*dy1+diff(x2j,y2)*dy2;
dx3j = diff(x3j,y1)*dy1+diff(x3j,y2)*dy2;
dx4j = diff(x4j,y1)*dy1+diff(x4j,y2)*dy2;
% Dissociation Constant Initialization Parameters
x2jk = (x4j^2*P)/(K2^2);
x3jk = (x1j^2)/(x2j^2*K1^2*P);
%solving the combined equations for dy1 and dy2
eqn = [diff(x2j,y1)*dy1+diff(x2j,y2)*dy2 == diff(x2jk,y1)*dy1 + diff(x2jk,y2)*dy2 + diff(x2jk,P)*dP + diff(x2jk,K1)*dK1+diff(x2jk,K2)*dK2, diff(x3j,y1)*dy1+diff(x3j,y2)*dy2 == diff(x3jk,y1)*dy1 + diff(x3jk,y2)*dy2 + diff(x3jk,P)*dP + diff(x3jk,K1)*dK1+diff(x3jk,K2)*dK2];
sol = solve(eqn,[dy1 dy2]);
dy1 = simplify(sol.dy1,"Steps",50);
dy2 = simplify(sol.dy2,"Steps",50);
%Now that we have equations for dy1 and dy2, we can proceed with the
%numerical analysis of the system
clearvars
h = 0.01;
tf = 60;
N = tf/h;
t = linspace(0,tf,N+1);
y1 = zeros(1,N+1);
y2 = zeros(1,N+1);
y1(1) = 0;
y2(1) = 0;
dK1 = 1;
dK2 = 1;
K1 = 1;
K2 = 1;
dP = 1;
P = 1;
for n=1:N
k1y1 = h*f(y1(n),K1,K2,dK1,dK2,dP,P,y1,y2);
k2y1 = h*f(y1(n)+1/2*k1y1,K1,K2,dK1,dK2,dP,P,y1,y2);
k3y1 = h*f(y1(n)+1/2*k2y1,K1,K2,dK1,dK2,dP,P,y1,y2);
k4y1 = h*f(y1(n)+k3y1,K1,K2,dK1,dK2,dP,P,y1,y2);
y1(n+1) = y1(n)+(1/6)*(k1y1+2*(k2y1+ k3y1)+ k4y1);
end
for n=1:N
k1y2 = h*g(y2(n),K1,K2,dK1,dK2,dP,P,y1,y2);
k2y2 = h*g(y2(n)+1/2*k1y2,K1,K2,dK1,dK2,dP,P,y1,y2);
k3y2 = h*g(y2(n)+1/2*k2y2,K1,K2,dK1,dK2,dP,P,y1,y2);
k4y2 = h*g(y2(n)+k3y2,K1,K2,dK1,dK2,dP,P,y1,y2);
y2(n+1) = y2(n)+(1/6)*(k1y2+2*(k2y2+k3y2)+k4y2);
end
nargin('dy1')
function dy1 = f(K1,K2,dK1,dK2,dP,P,y1,y2)
dy1 = (dP*K1^3*K2*P^2*y1^5*y2^4 - 4*dP*K1^3*K2*P^2*y1^5*y2^3 + 6*dP*K1^3*K2*P^2*y1^5*y2^2 - 4*dP*K1^3*K2*P^2*y1^5*y2 + dP*K1^3*K2*P^2*y1^5 + dP*K1^3*K2*P^2*y1^4*y2^4 - 4*dP*K1^3*K2*P^2*y1^4*y2^3 + 6*dP*K1^3*K2*P^2*y1^4*y2^2 - 4*dP*K1^3*K2*P^2*y1^4*y2 + dP*K1^3*K2*P^2*y1^4 - 2*dK2*K1^3*P^3*y1^5*y2^4 + 8*dK2*K1^3*P^3*y1^5*y2^3 - 12*dK2*K1^3*P^3*y1^5*y2^2 + 8*dK2*K1^3*P^3*y1^5*y2 - 2*dK2*K1^3*P^3*y1^5 - 2*dK2*K1^3*P^3*y1^4*y2^4 + 8*dK2*K1^3*P^3*y1^4*y2^3 - 12*dK2*K1^3*P^3*y1^4*y2^2 + 8*dK2*K1^3*P^3*y1^4*y2 - 2*dK2*K1^3*P^3*y1^4 + 4*dP*K1*K2^3*y1^5*y2^2 + 8*dP*K1*K2^3*y1^5*y2 + 4*dP*K1*K2^3*y1^5 - 4*dP*K1*K2^3*y1^4*y2^2 + 8*dP*K1*K2^3*y1^4*y2 + 12*dP*K1*K2^3*y1^4 - 4*dP*K1*K2^3*y1^3*y2^2 - 24*dP*K1*K2^3*y1^3*y2 - 4*dP*K1*K2^3*y1^3 + 4*dP*K1*K2^3*y1^2*y2^2 - 8*dP*K1*K2^3*y1^2*y2 - 28*dP*K1*K2^3*y1^2 + 16*dP*K1*K2^3*y1*y2 + 16*dP*K1*K2^3 - 4*dP*K1*K2*P*y1^5*y2^3 + 4*dP*K1*K2*P*y1^5*y2^2 + 4*dP*K1*K2*P*y1^5*y2 - 4*dP*K1*K2*P*y1^5 + 8*dP*K1*K2*P*y1^4*y2^3 - 16*dP*K1*K2*P*y1^4*y2^2 + 8*dP*K1*K2*P*y1^4*y2 - 4*dP*K1*K2*P*y1^3*y2^3 + 20*dP*K1*K2*P*y1^3*y2^2 - 28*dP*K1*K2*P*y1^3*y2 + 12*dP*K1*K2*P*y1^3 - 8*dP*K1*K2*P*y1^2*y2^2 + 16*dP*K1*K2*P*y1^2*y2 - 8*dP*K1*K2*P*y1^2 + 8*dK2*K1*P^2*y1^5*y2^3 + 8*dK2*K1*P^2*y1^5*y2^2 - 8*dK2*K1*P^2*y1^5*y2 - 8*dK2*K1*P^2*y1^5 - 16*dK2*K1*P^2*y1^4*y2^3 + 16*dK2*K1*P^2*y1^4*y2^2 + 16*dK2*K1*P^2*y1^4*y2 - 16*dK2*K1*P^2*y1^4 + 8*dK2*K1*P^2*y1^3*y2^3 - 56*dK2*K1*P^2*y1^3*y2^2 + 24*dK2*K1*P^2*y1^3*y2 + 24*dK2*K1*P^2*y1^3 + 32*dK2*K1*P^2*y1^2*y2^2 - 64*dK2*K1*P^2*y1^2*y2 + 32*dK2*K1*P^2*y1^2 + 32*dK2*K1*P^2*y1*y2 - 32*dK2*K1*P^2*y1 + 8*dK1*K2^3*P*y1^5*y2^2 + 16*dK1*K2^3*P*y1^5*y2 + 8*dK1*K2^3*P*y1^5 - 8*dK1*K2^3*P*y1^4*y2^2 + 16*dK1*K2^3*P*y1^4*y2 + 24*dK1*K2^3*P*y1^4 - 8*dK1*K2^3*P*y1^3*y2^2 - 48*dK1*K2^3*P*y1^3*y2 - 8*dK1*K2^3*P*y1^3 + 8*dK1*K2^3*P*y1^2*y2^2 - 16*dK1*K2^3*P*y1^2*y2 - 56*dK1*K2^3*P*y1^2 + 32*dK1*K2^3*P*y1*y2 + 32*dK1*K2^3*P + 16*dK1*K2*P^2*y1^5*y2^2 - 16*dK1*K2*P^2*y1^5 - 16*dK1*K2*P^2*y1^4*y2^2 + 32*dK1*K2*P^2*y1^4*y2 - 16*dK1*K2*P^2*y1^4 - 16*dK1*K2*P^2*y1^3*y2^2 - 32*dK1*K2*P^2*y1^3*y2 + 48*dK1*K2*P^2*y1^3 + 16*dK1*K2*P^2*y1^2*y2^2 - 32*dK1*K2*P^2*y1^2*y2 + 16*dK1*K2*P^2*y1^2 + 32*dK1*K2*P^2*y1*y2 - 32*dK1*K2*P^2*y1)/(16*K1*K2*P*(y1 - 1)*(y1 + y1*y2 + 2)*(K2^2*y1 - 2*P*y1 + 2*K2^2 + 2*P*y1*y2 + K2^2*y1*y2));
end
function dy2 = g(K1,K2,dK1,dK2,dP,P,y1,y2)
dy2 = -((y2 - 1)*(dP*K1^3*K2*P^2*y1^4*y2^4 - 4*dP*K1^3*K2*P^2*y1^4*y2^3 + 6*dP*K1^3*K2*P^2*y1^4*y2^2 - 4*dP*K1^3*K2*P^2*y1^4*y2 + dP*K1^3*K2*P^2*y1^4 - 2*dK2*K1^3*P^3*y1^4*y2^4 + 8*dK2*K1^3*P^3*y1^4*y2^3 - 12*dK2*K1^3*P^3*y1^4*y2^2 + 8*dK2*K1^3*P^3*y1^4*y2 - 2*dK2*K1^3*P^3*y1^4 + 4*dP*K1*K2^3*y1^4*y2^2 + 8*dP*K1*K2^3*y1^4*y2 + 4*dP*K1*K2^3*y1^4 - 8*dP*K1*K2^3*y1^3*y2^2 + 8*dP*K1*K2^3*y1^3 + 4*dP*K1*K2^3*y1^2*y2^2 - 24*dP*K1*K2^3*y1^2*y2 - 12*dP*K1*K2^3*y1^2 + 16*dP*K1*K2^3*y1*y2 - 16*dP*K1*K2^3*y1 + 16*dP*K1*K2^3 + 4*dP*K1*K2*P*y1^4*y2^3 + 12*dP*K1*K2*P*y1^4*y2^2 - 4*dP*K1*K2*P*y1^4*y2 - 12*dP*K1*K2*P*y1^4 - 4*dP*K1*K2*P*y1^3*y2^3 - 4*dP*K1*K2*P*y1^3*y2^2 + 20*dP*K1*K2*P*y1^3*y2 - 12*dP*K1*K2*P*y1^3 - 8*dP*K1*K2*P*y1^2*y2^2 - 16*dP*K1*K2*P*y1^2*y2 + 24*dP*K1*K2*P*y1^2 - 8*dK2*K1*P^2*y1^4*y2^3 - 8*dK2*K1*P^2*y1^4*y2^2 + 8*dK2*K1*P^2*y1^4*y2 + 8*dK2*K1*P^2*y1^4 + 8*dK2*K1*P^2*y1^3*y2^3 - 24*dK2*K1*P^2*y1^3*y2^2 - 8*dK2*K1*P^2*y1^3*y2 + 24*dK2*K1*P^2*y1^3 + 32*dK2*K1*P^2*y1^2*y2^2 - 32*dK2*K1*P^2*y1^2*y2 + 32*dK2*K1*P^2*y1*y2 - 32*dK2*K1*P^2*y1 + 8*dK1*K2^3*P*y1^4*y2^2 + 16*dK1*K2^3*P*y1^4*y2 + 8*dK1*K2^3*P*y1^4 - 16*dK1*K2^3*P*y1^3*y2^2 + 16*dK1*K2^3*P*y1^3 + 8*dK1*K2^3*P*y1^2*y2^2 - 48*dK1*K2^3*P*y1^2*y2 - 24*dK1*K2^3*P*y1^2 + 32*dK1*K2^3*P*y1*y2 - 32*dK1*K2^3*P*y1 + 32*dK1*K2^3*P + 16*dK1*K2*P^2*y1^4*y2^2 - 16*dK1*K2*P^2*y1^4 - 32*dK1*K2*P^2*y1^3*y2^2 + 32*dK1*K2*P^2*y1^3*y2 + 16*dK1*K2*P^2*y1^2*y2^2 - 64*dK1*K2*P^2*y1^2*y2 + 48*dK1*K2*P^2*y1^2 + 32*dK1*K2*P^2*y1*y2 - 32*dK1*K2*P^2*y1))/(16*K1*K2*P*y1*(y1 - 1)*(y1 + y1*y2 + 2)*(K2^2*y1 - 2*P*y1 + 2*K2^2 + 2*P*y1*y2 + K2^2*y1*y2));
end
13 Comments
@Jason Louison: "where the command window gives a "Too many input arguments" error for the function of dy1" - Do not improve a failing function, but fix the error at first. Whenever you mention an error in the forum, attach a copy of the complete message. Then the readers do not have to guess the details.
For the vectorization, replacing all * by .* and ^ by .^ should be sufficient already.
Jason Louison
on 20 Dec 2022
Torsten
on 20 Dec 2022
You can't solve for y1 and y2 separately - the differential equations are not independent from each other.
Jason Louison
on 20 Dec 2022
Edited: Jason Louison
on 20 Dec 2022
I am not solving them explicitly. I am solving the differential equations of their time derivates. By numerically integrating with the 4th Order Runge-Kutta method, I can use the rate of change of y1 and y2 to estimate the actual value of y1 and/or y2 at the next time step. I have done this with matlab before, and will attach a working code to demostrate:
% Chemical Equilibrium Reaction Solving: By Jason Louison
% Solve the following chemical equation: aO2 -> n1[O2]+n2[O]
% Goal: We want to be able to solve the chemical reaction above for all
% products assuming chemical equilibrium
% Assume a fraction of O2 dissociates. We will call this parameter "y"
%Write out the mass balance: a -> a(1-y)[O2]+2ay*phi[O]
%from this, we can derive the total number of products:
%np a(1-y)+2ay = a(1+y)
%From this, we define each mole fraction
% x1 = (1-y)/(1+y)
% x2 = (2y)/(1+y)
%From here, we can plug these values into the equilibrium constant
%equation and solve for the time derivate of y. To solve for the
%differential equation, initialization parameters xij are used.
%clearvars
syms K dK y dy P dP
x1i = (1-y)/(1+y);
x2i = (2*y)/(1+y);
dx1i = diff(x1i,y)*dy;
dx2i = diff(x2i,y)*dy;
x2iK = (x1i^2*P)/(K^2);
eqn = diff(x2i,y)*dy == diff(x2iK,y)*dy + diff(x2iK,P)*dP + diff(x2iK,K)*dK;
sol = solve(eqn,dy);
dy = sol;
%Simulation Parameters: Time step, final time, and Number of Cells
h = 0.01;
tf = 60;
N = tf/h;
%Designation of vector space for time and y
t = linspace(0,tf,N+1);
y = zeros(1,N+1);
%Initial Condition for y
y(1) = 0;
%Designation of vector space for dK, dP, P, and K:
dK = zeros(1,N+1);
dP = zeros(1,N+1);
P = zeros(1,N+1);
K = zeros(1,N+1);
%This allows us to vary the above parameters as functions of time. The
%parameters are assigned constant values below
dK = 0.34;
dP = 1;
P = 1;
K = 1;
for n=1:N %fourth-order runge kutta method algorithm
k1 = h*f(y(n),dK,dP,P,K);
k2 = h*f(y(n)+1/2*k1,dK,dP,P,K);
k3 = h*f(y(n)+1/2*k2,dK,dP,P,K);
k4 = h*f(y(n)+k3,dK,dP,P,K);
y(n+1) = y(n)+(1/6)*(k1+2*(k2+k3)+k4);
end
%Mole Fractions x1 and x2 as functions of dissociation fraction
x1 = (1-y)./(1+y);
x2 = (2.*y)./(1+y);
%Plot x1, x2, and y vs time
figure(1)
plot(t,x1,"r",t,x2,"b")
title('Dissociation of O_2 vs Time')
xlabel('Time (s)')
ylabel('mole fraction (ni/nt)')
legend({'O_2','O'})
%figure(2)
%plot(t,y)
%title('Dissociation Fraction of O_2 vs Time')
%xlabel('Time (s)')
%ylabel('mole fraction (ni/nt)')
%Main Differential Equation of y: This describes the rate of change of y
%with respect to itself, the time derivatives of Dissociation constant K
%and pressure P (dK and dP), as well as the dissociation constant K and
%pressure P:
function dy = f(y,dK,P,dP,K)
dy = (K.*dP - 2.*P.*dK - K.*dP*y + 2.*P.*dK.*y - K.*dP.*y.^2 + K.*dP.*y.^3 + 2.*P.*dK.*y.^2 - 2.*P.*dK.*y.^3)/(2.*K.^3.*y + 2.*K.^3 + 4.*K.*P - 4.*K.*P.*y);
end
Torsten
on 21 Dec 2022
Yes, but for the two-equation case, you might have noticed that dy1 depends on y2 and dy2 depends on y1. Thus the two equations cannot be solved independently from each other.
Jason Louison
on 21 Dec 2022
Yes. Imitate the user interface of the MATLAB ode solvers.
function dy = fun(t,y,pars)
Here, y = [y(1);y(2)] , time t and aybe additional parameters pars are inputs and [dy(1);dy(2)] are outputs.
Your Runge-Kutta code should address this functional interface. Thus in your updates
k1 = h*f(y(n),dK,dP,P,K);
k2 = h*f(y(n)+1/2*k1,dK,dP,P,K);
k3 = h*f(y(n)+1/2*k2,dK,dP,P,K);
k4 = h*f(y(n)+k3,dK,dP,P,K);
f,y(:,n) and ki will all be vectors.
Jason Louison
on 21 Dec 2022
Jason Louison
on 21 Dec 2022
Jan
on 21 Dec 2022
@Jason Louison: A general hint: Whenever you mention an error in the forum, post a copy of the complete message. Then the readers do not have to guess, which line is failing.
"Is there a way to do this in one step or command?" - seriously? The Replace dialog in the editor should do this sufficiently already. Even for really large functions with tenthousands of equations this takes fractions of a second only.
Using a fixed step Runge Kutta solver of order 4 solve a chemical reaction is a very bad idea. If the equation is stiff, the final value is dominated by the discretization error and more or less random. To reduce the noise, you have to run it with a tiny step size, such that the runtime will explode and the accumulated rounding errors will be the main part of the final value, such that you should not call it "result" anymore.
If the equation is stiff, you need a specific solver like ODE15s. Otherwise a solver with step-size control is the way to go, e.g. ODE45.
Solving an equation for a set of different parameters should not be done by vectorizing the function to be integrated, but by starting the integration in a loop, in which the parameters are varied. If multiple trajectories are calculated in the same integration, the component with the highest discretization and rounding errors determine the step-size. This means, that all other components are calculated too frequently and this reduces the accuracy due to the accumulation of rounding errors. The total run-time will increase because the system is integrated with a too small step-size in average.
If you have a really good reason to integrate a system with parameters provided as vectors, the variables to be integrated are not scalars anymore:
dy(1) = (dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^4 - 4.*dP.*K1.^3.*K2.*P.^2.*y(1).^5.*y(2).^3 + ...
dy(1) is a scalar, but the right hand side is a vector. Of course this cannot work. You have to blow up dy to a 2*(N+1) vector instead.
I still do not see the reason to vectorize the equations. This does not offer an advantage for integrating the system. The shown Runge-Kutta method is not sufficient to calculate the trajectories reliably in addition. This is not state of the art even in lessons for beginners, but an expensive random number generator with a poor entropy.
Use on of the built-in solvers with a step-size control and check the stiffness of the system (a cheap check: Try, if ODE45 or ODE15s is much faster).
Jason Louison
on 21 Dec 2022
Edited: Jason Louison
on 21 Dec 2022
@Jason Louison: Again: "the program kept giving me an “incorrect dimension error”." - If you mention an error, post a copy of the complete message and the relevant lines of code.
My comments are not harsh. I try to show you, how to solve the problem and how to avoid the typical errors made by beginners. Numerical maths is a demanding field of science. Feel free to ask questions about Matlab, because this is the purpose of this forum.
What is the actual problem you want to solve? You have a 2 DoF system of equations. Why do you want to solve the system "as an array"? What do you expect from vectorizing the function to be integrated?
One of your posted codes contain:
dK = zeros(1,N+1);
dP = zeros(1,N+1);
P = zeros(1,N+1);
K = zeros(1,N+1);
%This allows us to vary the above parameters as functions of time. The
%parameters are assigned constant values below
dK = 0.34;
dP = 1;
P = 1;
K = 1;
After creating dK etc. as vectors, you overwrite all of them by scalars. What is the underlying intention?
An example for a proper integration:
dK1 = 0;
dK2 = 0;
K1 = 0.5;
K2 = 0.7;
dP = 0;
P = 1;
[T, Y] = ode45(@(t,y) f(t,y, K1,K2,dK1,dK2,dP,P), [0, 60], [1, 0.5]);
plot(T, Y)
function dy = f(t,y, K1,K2,dK1,dK2,dP,P)
y1_5 = y(1)^5;
y1_4 = y(1)^4;
y1_3 = y(1)^3;
y1_2 = y(1)^2;
y2_4 = y(2)^4;
y2_3 = y(2)^3;
y2_2 = y(2)^2;
K2_3 = K2^3;
K2_2 = K2^2;
P_2 = P^2;
dy = [(dP*K1^3*K2*P_2*y1_5*y2_4 - 4*dP*K1^3*K2*P_2*y1_5*y2_3 + ...
6*dP*K1^3*K2*P_2*y1_5*y2_2 - 4*dP*K1^3*K2*P_2*y1_5*y(2) + ...
dP*K1^3*K2*P_2*y1_5 + dP*K1^3*K2*P_2*y1_4*y2_4 - ...
4*dP*K1^3*K2*P_2*y1_4*y2_3 + 6*dP*K1^3*K2*P_2*y1_4*y2_2 - ...
4*dP*K1^3*K2*P_2*y1_4*y(2) + dP*K1^3*K2*P_2*y1_4 - ...
2*dK2*K1^3*P^3*y1_5*y2_4 + 8*dK2*K1^3*P^3*y1_5*y2_3 - ...
12*dK2*K1^3*P^3*y1_5*y2_2 + 8*dK2*K1^3*P^3*y1_5*y(2) - ...
2*dK2*K1^3*P^3*y1_5 - 2*dK2*K1^3*P^3*y1_4*y2_4 + ...
8*dK2*K1^3*P^3*y1_4*y2_3 - 12*dK2*K1^3*P^3*y1_4*y2_2 + ...
8*dK2*K1^3*P^3*y1_4*y(2) - 2*dK2*K1^3*P^3*y1_4 + ...
4*dP*K1*K2_3*y1_5*y2_2 + 8*dP*K1*K2_3*y1_5*y(2) + ...
4*dP*K1*K2_3*y1_5 - 4*dP*K1*K2_3*y1_4*y2_2 + ...
8*dP*K1*K2_3*y1_4*y(2) + 12*dP*K1*K2_3*y1_4 - ...
4*dP*K1*K2_3*y1_3*y2_2 - 24*dP*K1*K2_3*y1_3*y(2) - ...
4*dP*K1*K2_3*y1_3 + 4*dP*K1*K2_3*y1_2*y2_2 - ...
8*dP*K1*K2_3*y1_2*y(2) - 28*dP*K1*K2_3*y1_2 + ...
16*dP*K1*K2_3*y(1)*y(2) + 16*dP*K1*K2_3 - ...
4*dP*K1*K2*P*y1_5*y2_3 + 4*dP*K1*K2*P*y1_5*y2_2 + ...
4*dP*K1*K2*P*y1_5*y(2) - 4*dP*K1*K2*P*y1_5 + ...
8*dP*K1*K2*P*y1_4*y2_3 - 16*dP*K1*K2*P*y1_4*y2_2 + ...
8*dP*K1*K2*P*y1_4*y(2) - 4*dP*K1*K2*P*y1_3*y2_3 + ...
20*dP*K1*K2*P*y1_3*y2_2 - 28*dP*K1*K2*P*y1_3*y(2) + ...
12*dP*K1*K2*P*y1_3 - 8*dP*K1*K2*P*y1_2*y2_2 + ...
16*dP*K1*K2*P*y1_2*y(2) - 8*dP*K1*K2*P*y1_2 + ...
8*dK2*K1*P_2*y1_5*y2_3 + 8*dK2*K1*P_2*y1_5*y2_2 - ...
8*dK2*K1*P_2*y1_5*y(2) - 8*dK2*K1*P_2*y1_5 - ...
16*dK2*K1*P_2*y1_4*y2_3 + 16*dK2*K1*P_2*y1_4*y2_2 + ...
16*dK2*K1*P_2*y1_4*y(2) - 16*dK2*K1*P_2*y1_4 + ...
8*dK2*K1*P_2*y1_3*y2_3 - 56*dK2*K1*P_2*y1_3*y2_2 + ...
24*dK2*K1*P_2*y1_3*y(2) + 24*dK2*K1*P_2*y1_3 + ...
32*dK2*K1*P_2*y1_2*y2_2 - 64*dK2*K1*P_2*y1_2*y(2) + ...
32*dK2*K1*P_2*y1_2 + 32*dK2*K1*P_2*y(1)*y(2) - ...
32*dK2*K1*P_2*y(1) + 8*dK1*K2_3*P*y1_5*y2_2 + ...
16*dK1*K2_3*P*y1_5*y(2) + 8*dK1*K2_3*P*y1_5 - ...
8*dK1*K2_3*P*y1_4*y2_2 + 16*dK1*K2_3*P*y1_4*y(2) + ...
24*dK1*K2_3*P*y1_4 - 8*dK1*K2_3*P*y1_3*y2_2 - ...
48*dK1*K2_3*P*y1_3*y(2) - 8*dK1*K2_3*P*y1_3 + ...
8*dK1*K2_3*P*y1_2*y2_2 - 16*dK1*K2_3*P*y1_2*y(2) - ...
56*dK1*K2_3*P*y1_2 + 32*dK1*K2_3*P*y(1)*y(2) + ...
32*dK1*K2_3*P + 16*dK1*K2*P_2*y1_5*y2_2 - ...
16*dK1*K2*P_2*y1_5 - 16*dK1*K2*P_2*y1_4*y2_2 + ...
32*dK1*K2*P_2*y1_4*y(2) - 16*dK1*K2*P_2*y1_4 - ...
16*dK1*K2*P_2*y1_3*y2_2 - 32*dK1*K2*P_2*y1_3*y(2) + ...
48*dK1*K2*P_2*y1_3 + 16*dK1*K2*P_2*y1_2*y2_2 - ...
32*dK1*K2*P_2*y1_2*y(2) + 16*dK1*K2*P_2*y1_2 + ...
32*dK1*K2*P_2*y(1)*y(2) - 32*dK1*K2*P_2*y(1)) / ...
(16*K1*K2*P*(y(1) - 1)*(y(1) + y(1)*y(2) + 2)*(K2_2*y(1) - ...
2*P*y(1) + 2*K2_2 + 2*P*y(1)*y(2) + K2_2*y(1)*y(2))); ...
-((y(2) - 1)*(dP*K1^3*K2*P_2*y1_4*y2_4 - 4*dP*K1^3*K2*P_2*y1_4*y2_3 + ...
6*dP*K1^3*K2*P_2*y1_4*y2_2 - 4*dP*K1^3*K2*P_2*y1_4*y(2) + ...
dP*K1^3*K2*P_2*y1_4 - 2*dK2*K1^3*P^3*y1_4*y2_4 + ...
8*dK2*K1^3*P^3*y1_4*y2_3 - 12*dK2*K1^3*P^3*y1_4*y2_2 + ...
8*dK2*K1^3*P^3*y1_4*y(2) - 2*dK2*K1^3*P^3*y1_4 + ...
4*dP*K1*K2_3*y1_4*y2_2 + 8*dP*K1*K2_3*y1_4*y(2) + ...
4*dP*K1*K2_3*y1_4 - 8*dP*K1*K2_3*y1_3*y2_2 + ...
8*dP*K1*K2_3*y1_3 + 4*dP*K1*K2_3*y1_2*y2_2 - ...
24*dP*K1*K2_3*y1_2*y(2) - 12*dP*K1*K2_3*y1_2 + ...
16*dP*K1*K2_3*y(1)*y(2) - 16*dP*K1*K2_3*y(1) + ...
16*dP*K1*K2_3 + 4*dP*K1*K2*P*y1_4*y2_3 + ...
12*dP*K1*K2*P*y1_4*y2_2 - 4*dP*K1*K2*P*y1_4*y(2) - ...
12*dP*K1*K2*P*y1_4 - 4*dP*K1*K2*P*y1_3*y2_3 - ...
4*dP*K1*K2*P*y1_3*y2_2 + 20*dP*K1*K2*P*y1_3*y(2) - ...
12*dP*K1*K2*P*y1_3 - 8*dP*K1*K2*P*y1_2*y2_2 - ...
16*dP*K1*K2*P*y1_2*y(2) + 24*dP*K1*K2*P*y1_2 - ...
8*dK2*K1*P_2*y1_4*y2_3 - 8*dK2*K1*P_2*y1_4*y2_2 + ...
8*dK2*K1*P_2*y1_4*y(2) + 8*dK2*K1*P_2*y1_4 + ...
8*dK2*K1*P_2*y1_3*y2_3 - 24*dK2*K1*P_2*y1_3*y2_2 - ...
8*dK2*K1*P_2*y1_3*y(2) + 24*dK2*K1*P_2*y1_3 + ...
32*dK2*K1*P_2*y1_2*y2_2 - 32*dK2*K1*P_2*y1_2*y(2) + ...
32*dK2*K1*P_2*y(1)*y(2) - 32*dK2*K1*P_2*y(1) + ...
8*dK1*K2_3*P*y1_4*y2_2 + 16*dK1*K2_3*P*y1_4*y(2) + ...
8*dK1*K2_3*P*y1_4 - 16*dK1*K2_3*P*y1_3*y2_2 + ...
16*dK1*K2_3*P*y1_3 + 8*dK1*K2_3*P*y1_2*y2_2 - ...
48*dK1*K2_3*P*y1_2*y(2) - 24*dK1*K2_3*P*y1_2 + ...
32*dK1*K2_3*P*y(1)*y(2) - 32*dK1*K2_3*P*y(1) + ...
32*dK1*K2_3*P + 16*dK1*K2*P_2*y1_4*y2_2 - ...
16*dK1*K2*P_2*y1_4 - 32*dK1*K2*P_2*y1_3*y2_2 + ...
32*dK1*K2*P_2*y1_3*y(2) + 16*dK1*K2*P_2*y1_2*y2_2 - ...
64*dK1*K2*P_2*y1_2*y(2) + 48*dK1*K2*P_2*y1_2 + ...
32*dK1*K2*P_2*y(1)*y(2) - 32*dK1*K2*P_2*y(1))) / ...
(16*K1*K2*P*y(1)*(y(1) - 1)*(y(1) + y(1)*y(2) + 2)*(K2_2*y(1) - ...
2*P*y(1) + 2*K2_2 + 2*P*y(1)*y(2) + K2_2*y(1)*y(2)))];
end
The output is empty, because the calculated trajectory contains NaN values only. I assume it is divided by zero.
So a vectorization is not meaningful at all. It is not clear yet, what you want to achieve and the posted functions does not have a useful result.
Jason Louison
on 21 Dec 2022
Accepted Answer
See how you can use the below code for your purpose.
The Runge-Kutta code is made to deal with systems of differential equations of arbitrary size (thus with vector inputs).
tstart = 0.0;
tend = 1.0;
nt = 100;
h = (tend - tstart)/nt;
T = (tstart:h:tend).';
Y0 = [1 -1];
B = 4;
Y = runge_kutta_RK4(@(t,y)f(t,y,B),T,Y0);
plot(T,Y)
function Y = runge_kutta_RK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for i = 2:N
t = T(i-1);
y = Y(i-1,:);
h = T(i) - T(i-1);
k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
end
function dy = f(t,y,B)
dy = zeros(1,2);
dy(1) = y(2);
dy(2) = -exp(-B*t)-y(1)+5*exp(-2*t)-2*exp(-(B+2)*t)+exp(-B*t)+t;
end
13 Comments
Jason Louison
on 21 Dec 2022
Jason Louison
on 22 Dec 2022
Edited: Jason Louison
on 22 Dec 2022
Tortsen,
After hours of trial and error and research, I have finally got the program to work. It does exactly what I want it to do. I cannot thank you and Jan enough. The latest version of the code is shown below. Attached is the output for the given input parameters and initial values.
I do believe the program can be optimized even more, especially with the use of the ode45 solver, but I wanted to see a runge kutta 4th order algorithm be put to use. Other areas of improvement could be a more organized layout as well as a less confusing variable input system (The functions MUST have their own variables, and the values I assign within the function don't get used. I find it strange). Another worry of mine is that the initial conditions cannot be zero for both y(1) and y(2). One has to be greater than zero, and calculating this right before the reaction begins may prove difficult. The main flaw of the program that was keeping me from being able to solve the system was that I needed a for loop to tell the differential equations what the value of each parameter was each time step. Nevertheless, I am extremely happy and grateful, as this code serves as a basis for solving more complex reactions. Here is a result figure I obtained from a different run, now showing the amount of each substance as well as dissociation and mole fraction:
% Chemical Equilibrium
% Solve the following chemical equation:
% aCO2 -> n1[CO2]+n2[O2]+n3[CO]+n4[O]
% Goal: We want to be able to solve the chemical reaction above for all
% products assuming chemical equilibrium
% Assume a fraction of CO2 and O2 dissociates. We will call these parameters
% "y1" and "y2"
% Write out the mass balance:
% a -> a(1-y1)[CO2]+(1/2)ay1(1-y2)[O2]+ay1[CO]+ay1[CO]+ay1y2[O]
% From this, we can derive the total number of products:
% np = a[1-y1+(1/2)y1(1-y2)+y1+y1y2]
% np = (1/2)a[2+y1(1+y2)]
% From this, we define each mole fraction
% x1 = 2(1-y1)/(2+y1(1+y2))
% x2 = y1(1-y2)/(2+y1(1+y2))
% x3 = 2y1/(2+y1(1+y2))
% x4 = 2y1y2/(2+y1(1+y2))
% From here, we can plug these values into the equilibrium constant
% equation and solve for the time derivate of y. To solve for the
% differential equation, initialization parameters xij are used.
%clearvars
% Initialization Parameters: These mole fraction and dissociation constant
% equations are being used to solve for the time derivatives of y1 and y2
clearvars
syms K1 dK1 K2 dK2 y1 dy1 y2 dy2 P dP
x1j = 2*(1-y1)/(2+y1*(1+y2));
x2j = y1*(1-y2)/(2+y1*(1+y2));
x3j = 2*y1/(2+y1*(1+y2));
x4j = (2*y1*y2)/(2+y1*(1+y2));
% Differential Initialization Parameters
dx1j = diff(x1j,y1)*dy1+diff(x1j,y2)*dy2;
dx2j = diff(x2j,y1)*dy1+diff(x2j,y2)*dy2;
dx3j = diff(x3j,y1)*dy1+diff(x3j,y2)*dy2;
dx4j = diff(x4j,y1)*dy1+diff(x4j,y2)*dy2;
% Dissociation Constant Initialization Parameters
x2jk = (x4j^2*P)/(K2^2);
x3jk = (x1j^2)/(x2j^2*K1^2*P);
%solving the combined equations for dy1 and dy2
eqn = [(diff(x2j,y1)*dy1 + diff(x2j,y2)*dy2 == diff(x2jk,y1)*dy1 + diff(x2jk,y2)*dy2 + diff(x2jk,P)*dP + diff(x2jk,K1)*dK1 + diff(x2jk,K2)*dK2)
(diff(x3j,y1)*dy1 + diff(x3j,y2)*dy2 == diff(x3jk,y1)*dy1 + diff(x3jk,y2)*dy2 + diff(x3jk,P)*dP + diff(x3jk,K1)*dK1 + diff(x3jk,K2)*dK2)];
sol = solve(eqn,[dy1 dy2]);
dy1i = simplify(sol.dy1,"Steps",50);
dy2i = simplify(sol.dy2,"Steps",50);
%Copy and Paste the equations from the variables dy1i and dy2i into the
%differential equation functions at the bottom of the code. Be sure to
%replace all operators with array compatible ones, all parametes with a t
%at the end, and y1 and y2 to y(1) and y(2)
%Now that we have equations for dy1 and dy2, we can proceed with the
%numerical analysis of the system
clearvars
tstart = 0.0;
tend = 120.0;
nt = 1000000;
h = (tend - tstart)/nt;
T = (tstart:h:tend).';
%Values of Parameters. These are stored in the matrix elements for each. It
%can be a function of time or a constant value.
%Global Variables
dK1 = sin(3.14*T);
dK2 = 2*cos(3.14*T);
K1 = 60;
K2 = 25;
dP = 1;
P = 1;
%Initial Conditions for Dissociation fraction y1 and y2
Y0 = [0.01 0.00];
%Array variables defining global varables at time t as an array
dK1f = dK1 + zeros(nt+1,1);
dK2f = dK2 + zeros(nt+1,1);
K1f = K1 + zeros(nt+1,1);
K2f = K2 + zeros(nt+1,1);
dPf = dP + zeros(nt+1,1);
Pf = P + zeros(nt+1,1);
yf = zeros(nt+1,2);
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for i = 2:N
dK1t = dK1f(i-1);
dK2t = dK2f(i-1);
K1t = K1f(i-1);
K2t = K2f(i-1);
dPt = dPf(i-1);
Pt = Pf(i-1);
end
Y = runge_kutta_RK4(@(t,y)f(t,K1t,K2t,dK1t,dK2t,dPt,Pt,y),T,Y0);
x1 = 2.*(1-Y(:,1))./(2+Y(:,1).*(1+Y(:,2)));
x2 = (Y(:,1).*(1-Y(:,2)))./(2+Y(:,1).*(1+Y(:,2)));
x3 = 2.*Y(:,1)./(2+Y(:,1).*(1+Y(:,2)));
x4 = (2.*Y(:,1).*Y(:,2))./(2+Y(:,1).*(1+Y(:,2)));
subplot(1,2,1);
plot(T,Y(:,1),T,Y(:,2))
title('Dissociation Fraction of CO_2 and O_2 vs. Time')
xlabel('Time (s)')
ylabel('Dissociation Fraction (ni/nt)')
legend({'CO_2','O_2'})
subplot(1,2,2)
plot(T,x1,T,x2,T,x3,T,x4)
title('Mole Fractions of CO_2, O_2, CO, and O vs. Time')
xlabel('Time (s)')
ylabel('Mole Fraction (ni/nt)')
legend({'CO_2','O_2','CO','O'})
% Differential Equation Functions
function dy = f(~,dK1t,dK2t,K1t,K2t,dPt,Pt,y)
dy = zeros(1,2);
dy(1) = -(2.*(- dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^5 + 3.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^4 - 3.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^3 + ...
dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^2 + 2.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^5 - 6.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^4 + 6.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^3 - ...
2.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^5.*y(2).^3 - dPt.*K1t.*K2t.^3.*y(1).^5.*y(2).^2 + dPt.*K1t.*K2t.^3.*y(1).^5.*y(2) + dPt.*K1t.*K2t.^3.*y(1).^5 + ...
dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^3 - 3.*dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^4.*y(2) + 3.*dPt.*K1t.*K2t.^3.*y(1).^4 + dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^3 + ...
5.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^2 - 5.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2) - dPt.*K1t.*K2t.^3.*y(1).^3 - dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^3 + 3.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^2 + ...
5.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2) - 7.*dPt.*K1t.*K2t.^3.*y(1).^2 - 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2).^2 + 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2) - 4.*dPt.*K1t.*K2t.^3.*y(2) + ...
4.*dPt.*K1t.*K2t.^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^4 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^2 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2) - ...
8.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^4 + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2) + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^4 - 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^3 - ...
12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^2 - 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2) + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^3 + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^2 - ...
16.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2) + 16.*dPt.*K1t.*K2t.*Pt.*y(1).*y(2) - 8.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^3 - ...
8.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^2 + 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^2 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^4 + ...
48.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^3 + 24.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^2 - 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^3 + 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^2 - ...
32.*dK2t.*K1t.*Pt.^2.*y(1).*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2).^3 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2).^2 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2) + ...
2.*dK1t.*K2t.^3.*Pt.*y(1).^5 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^3 - 6.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2) + 6.*dK1t.*K2t.^3.*Pt.*y(1).^4 + ...
2.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^3 + 10.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^2 - 10.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2) - 2.*dK1t.*K2t.^3.*Pt.*y(1).^3 - ...
2.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^3 + 6.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^2 + 10.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2) - 14.*dK1t.*K2t.^3.*Pt.*y(1).^2 - ...
8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2).^2 + 8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2) - 8.*dK1t.*K2t.^3.*Pt.*y(2) + 8.*dK1t.*K2t.^3.*Pt - 8.*dK1t.*K2t.*Pt.^2.*y(1).^5.*y(2).^3 + ...
8.*dK1t.*K2t.*Pt.^2.*y(1).^5.*y(2) - 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2) + 24.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2).^3 - ...
24.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2) - 16.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^3 + 48.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^2 - 32.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2) - ...
32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2).^2 + 32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2)))/(K1t.*K2t.*Pt.*(- K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^4 + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^3 - ...
2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2) + K1t.^2.*K2t.^2.*Pt.*y(1).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^3 + 6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^2 - ...
6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2) + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2 - 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^4 + 24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^3 - ...
24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^2 + 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2) + 8.*K2t.^2.*y(1).^3.*y(2).^3 + 8.*K2t.^2.*y(1).^3.*y(2).^2 - 8.*K2t.^2.*y(1).^3.*y(2) - ...
8.*K2t.^2.*y(1).^3 - 8.*K2t.^2.*y(1).^2.*y(2).^3 + 24.*K2t.^2.*y(1).^2.*y(2).^2 + 8.*K2t.^2.*y(1).^2.*y(2) - 24.*K2t.^2.*y(1).^2 - ...
32.*K2t.^2.*y(1).*y(2).^2 + 32.*K2t.^2.*y(1).*y(2) - 32.*K2t.^2.*y(2) + 32.*K2t.^2 + 48.*Pt.*y(1).^3.*y(2).^3 + 32.*Pt.*y(1).^3.*y(2).^2 - 16.*Pt.*y(1).^3.*y(2) - ...
48.*Pt.*y(1).^2.*y(2).^3 + 32.*Pt.*y(1).^2.*y(2).^2 - 48.*Pt.*y(1).^2.*y(2) - 64.*Pt.*y(1).*y(2).^2 + 64.*Pt.*y(2)));
dy(2) = -(2.*(y(2) - 1).*(- 2.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^4 + 4.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^3 - 2.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 + ...
4.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^4 - 8.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^3 + 4.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^2 + dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^3 + ...
dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^4.*y(2) - dPt.*K1t.*K2t.^3.*y(1).^4 - 2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^3 + 2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^2 + ...
2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2) - 2.*dPt.*K1t.*K2t.^3.*y(1).^3 + dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^3 - 7.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^2 + ...
3.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2) + 3.*dPt.*K1t.*K2t.^3.*y(1).^2 + 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2).^2 - 8.*dPt.*K1t.*K2t.^3.*y(1).*y(2) + 4.*dPt.*K1t.*K2t.^3.*y(1) + ...
4.*dPt.*K1t.*K2t.^3.*y(2) - 4.*dPt.*K1t.*K2t.^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^4 + 16.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^3 + 12.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^2 - ...
4.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^4 - 8.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^3 + 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^2 - 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^3 - ...
24.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^2 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^3 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^2 + ...
8.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^3 - 24.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^2 + 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^3 + ...
32.*dK2t.*K1t.*Pt.^2.*y(1).*y(2).^2 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^3 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2) - ...
2.*dK1t.*K2t.^3.*Pt.*y(1).^4 - 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^3 + 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^2 + 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2) - ...
4.*dK1t.*K2t.^3.*Pt.*y(1).^3 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^3 - 14.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^2 + 6.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2) + ...
6.*dK1t.*K2t.^3.*Pt.*y(1).^2 + 8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2).^2 - 16.*dK1t.*K2t.^3.*Pt.*y(1).*y(2) + 8.*dK1t.*K2t.^3.*Pt.*y(1) + 8.*dK1t.*K2t.^3.*Pt.*y(2) - ...
8.*dK1t.*K2t.^3.*Pt + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^3 + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 - 32.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2).^3 + ...
16.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^3 - 48.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^2 + 32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2).^2))/(K1t.*K2t.*Pt.*y(1).*(- K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^4 + ...
2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2) + K1t.^2.*K2t.^2.*Pt.*y(1).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^3 + ...
6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^2 - 6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2) + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2 - 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^4 + ...
24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^3 - 24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^2 + 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2) + 8.*K2t.^2.*y(1).^3.*y(2).^3 + ...
8.*K2t.^2.*y(1).^3.*y(2).^2 - 8.*K2t.^2.*y(1).^3.*y(2) - 8.*K2t.^2.*y(1).^3 - 8.*K2t.^2.*y(1).^2.*y(2).^3 + 24.*K2t.^2.*y(1).^2.*y(2).^2 + ...
8.*K2t.^2.*y(1).^2.*y(2) - 24.*K2t.^2.*y(1).^2 - 32.*K2t.^2.*y(1).*y(2).^2 + 32.*K2t.^2.*y(1).*y(2) - 32.*K2t.^2.*y(2) + 32.*K2t.^2 + 48.*Pt.*y(1).^3.*y(2).^3 + ...
32.*Pt.*y(1).^3.*y(2).^2 - 16.*Pt.*y(1).^3.*y(2) - 48.*Pt.*y(1).^2.*y(2).^3 + 32.*Pt.*y(1).^2.*y(2).^2 - 48.*Pt.*y(1).^2.*y(2) - 64.*Pt.*y(1).*y(2).^2 + 64.*Pt.*y(2)));
dy = [dy(1); dy(2)];
end
function Y = runge_kutta_RK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
dK1f = 1;
dK2f = 1;
K1f = 1;
K2f = 1;
dPf = 1;
Pf = 1;
nt = 10000000;
dK1 = dK1f + zeros(nt+1,1);
dK2 = dK2f + zeros(nt+1,1);
K1 = K1f + zeros(nt+1,1);
K2 = K2f + zeros(nt+1,1);
dP = dPf + zeros(nt+1,1);
P = Pf + zeros(nt+1,1);
y = zeros(nt+1,2);
for i = 2:N
dK1t = dK1(i-1);
dK2t = dK2(i-1);
K1t = K1(i-1);
K2t = K2(i-1);
dPt = dP(i-1);
Pt = P(i-1);
t = T(i-1);
y = Y(i-1,:);
h = T(i) - T(i-1);
k0 = h*f(t,y)';
k1 = h*f(t,y+k0*0.5)';
k2 = h*f(t,y+k1*0.5)';
k3 = h*f(t,y+k2)';
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
end
I get different results.
I wonder why you modified the function "runge_kutta_RK4" and made it work incorrectly. Changes to "runge_kutta_RK4" are never necessary. Changes should to be made in the calling function and in the function where you supply the derivatives. Or do you modify MATLAB's inbuilt ODE45 when you want to solve a differential equation ?
tstart = 0.0;
tend = 120.0;
nt = 1000000;
h = (tend - tstart)/nt;
T = (tstart:h:tend).';
%Initial Conditions for Dissociation fraction y1 and y2
Y0 = [0.01 0.00];
N = numel(T);
n = numel(Y0);
Y = runge_kutta_RK4(@f,T,Y0);
x1 = 2.*(1-Y(:,1))./(2+Y(:,1).*(1+Y(:,2)));
x2 = (Y(:,1).*(1-Y(:,2)))./(2+Y(:,1).*(1+Y(:,2)));
x3 = 2.*Y(:,1)./(2+Y(:,1).*(1+Y(:,2)));
x4 = (2.*Y(:,1).*Y(:,2))./(2+Y(:,1).*(1+Y(:,2)));
subplot(1,2,1);
plot(T,Y(:,1),T,Y(:,2))
title('Dissociation Fraction of CO_2 and O_2 vs. Time')
xlabel('Time (s)')
ylabel('Dissociation Fraction (ni/nt)')
legend({'CO_2','O_2'})
subplot(1,2,2)
plot(T,x1,T,x2,T,x3,T,x4)
title('Mole Fractions of CO_2, O_2, CO, and O vs. Time')
xlabel('Time (s)')
ylabel('Mole Fraction (ni/nt)')
legend({'CO_2','O_2','CO','O'})
% Differential Equation Functions
function dy = f(t,y)
dK1t = sin(3.14*t);
dK2t = 2*cos(3.14*t);
K1t = 60;
K2t = 25;
dPt = 1;
Pt = 1;
dy = zeros(1,2);
dy(1) = -(2.*(- dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^5 + 3.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^4 - 3.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^3 + ...
dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^5.*y(2).^2 + 2.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^5 - 6.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^4 + 6.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^3 - ...
2.*dK2t.*K1t.^3.*Pt.^3.*y(1).^5.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^5.*y(2).^3 - dPt.*K1t.*K2t.^3.*y(1).^5.*y(2).^2 + dPt.*K1t.*K2t.^3.*y(1).^5.*y(2) + dPt.*K1t.*K2t.^3.*y(1).^5 + ...
dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^3 - 3.*dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^4.*y(2) + 3.*dPt.*K1t.*K2t.^3.*y(1).^4 + dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^3 + ...
5.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^2 - 5.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2) - dPt.*K1t.*K2t.^3.*y(1).^3 - dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^3 + 3.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^2 + ...
5.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2) - 7.*dPt.*K1t.*K2t.^3.*y(1).^2 - 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2).^2 + 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2) - 4.*dPt.*K1t.*K2t.^3.*y(2) + ...
4.*dPt.*K1t.*K2t.^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^4 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2).^2 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^5.*y(2) - ...
8.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^4 + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2) + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^4 - 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^3 - ...
12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^2 - 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2) + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^3 + 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^2 - ...
16.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2) + 16.*dPt.*K1t.*K2t.*Pt.*y(1).*y(2) - 8.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^3 - ...
8.*dK2t.*K1t.*Pt.^2.*y(1).^5.*y(2).^2 + 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^2 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^4 + ...
48.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^3 + 24.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^2 - 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^3 + 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^2 - ...
32.*dK2t.*K1t.*Pt.^2.*y(1).*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2).^3 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2).^2 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^5.*y(2) + ...
2.*dK1t.*K2t.^3.*Pt.*y(1).^5 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^3 - 6.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2) + 6.*dK1t.*K2t.^3.*Pt.*y(1).^4 + ...
2.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^3 + 10.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^2 - 10.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2) - 2.*dK1t.*K2t.^3.*Pt.*y(1).^3 - ...
2.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^3 + 6.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^2 + 10.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2) - 14.*dK1t.*K2t.^3.*Pt.*y(1).^2 - ...
8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2).^2 + 8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2) - 8.*dK1t.*K2t.^3.*Pt.*y(2) + 8.*dK1t.*K2t.^3.*Pt - 8.*dK1t.*K2t.*Pt.^2.*y(1).^5.*y(2).^3 + ...
8.*dK1t.*K2t.*Pt.^2.*y(1).^5.*y(2) - 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2) + 24.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2).^3 - ...
24.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2) - 16.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^3 + 48.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^2 - 32.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2) - ...
32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2).^2 + 32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2)))/(K1t.*K2t.*Pt.*(- K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^4 + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^3 - ...
2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2) + K1t.^2.*K2t.^2.*Pt.*y(1).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^3 + 6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^2 - ...
6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2) + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2 - 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^4 + 24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^3 - ...
24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^2 + 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2) + 8.*K2t.^2.*y(1).^3.*y(2).^3 + 8.*K2t.^2.*y(1).^3.*y(2).^2 - 8.*K2t.^2.*y(1).^3.*y(2) - ...
8.*K2t.^2.*y(1).^3 - 8.*K2t.^2.*y(1).^2.*y(2).^3 + 24.*K2t.^2.*y(1).^2.*y(2).^2 + 8.*K2t.^2.*y(1).^2.*y(2) - 24.*K2t.^2.*y(1).^2 - ...
32.*K2t.^2.*y(1).*y(2).^2 + 32.*K2t.^2.*y(1).*y(2) - 32.*K2t.^2.*y(2) + 32.*K2t.^2 + 48.*Pt.*y(1).^3.*y(2).^3 + 32.*Pt.*y(1).^3.*y(2).^2 - 16.*Pt.*y(1).^3.*y(2) - ...
48.*Pt.*y(1).^2.*y(2).^3 + 32.*Pt.*y(1).^2.*y(2).^2 - 48.*Pt.*y(1).^2.*y(2) - 64.*Pt.*y(1).*y(2).^2 + 64.*Pt.*y(2)));
dy(2) = -(2.*(y(2) - 1).*(- 2.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^4 + 4.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^3 - 2.*dPt.*K1t.^3.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 + ...
4.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^4 - 8.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^3 + 4.*dK2t.*K1t.^3.*Pt.^3.*y(1).^4.*y(2).^2 + dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^3 + ...
dPt.*K1t.*K2t.^3.*y(1).^4.*y(2).^2 - dPt.*K1t.*K2t.^3.*y(1).^4.*y(2) - dPt.*K1t.*K2t.^3.*y(1).^4 - 2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^3 + 2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2).^2 + ...
2.*dPt.*K1t.*K2t.^3.*y(1).^3.*y(2) - 2.*dPt.*K1t.*K2t.^3.*y(1).^3 + dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^3 - 7.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2).^2 + ...
3.*dPt.*K1t.*K2t.^3.*y(1).^2.*y(2) + 3.*dPt.*K1t.*K2t.^3.*y(1).^2 + 4.*dPt.*K1t.*K2t.^3.*y(1).*y(2).^2 - 8.*dPt.*K1t.*K2t.^3.*y(1).*y(2) + 4.*dPt.*K1t.*K2t.^3.*y(1) + ...
4.*dPt.*K1t.*K2t.^3.*y(2) - 4.*dPt.*K1t.*K2t.^3 + 4.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^4 + 16.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^3 + 12.*dPt.*K1t.*K2t.*Pt.*y(1).^4.*y(2).^2 - ...
4.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^4 - 8.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^3 + 12.*dPt.*K1t.*K2t.*Pt.*y(1).^3.*y(2).^2 - 8.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^3 - ...
24.*dPt.*K1t.*K2t.*Pt.*y(1).^2.*y(2).^2 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^3 - 8.*dK2t.*K1t.*Pt.^2.*y(1).^4.*y(2).^2 + ...
8.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^4 - 16.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^3 - 24.*dK2t.*K1t.*Pt.^2.*y(1).^3.*y(2).^2 + 32.*dK2t.*K1t.*Pt.^2.*y(1).^2.*y(2).^3 + ...
32.*dK2t.*K1t.*Pt.^2.*y(1).*y(2).^2 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^3 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2).^2 - 2.*dK1t.*K2t.^3.*Pt.*y(1).^4.*y(2) - ...
2.*dK1t.*K2t.^3.*Pt.*y(1).^4 - 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^3 + 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2).^2 + 4.*dK1t.*K2t.^3.*Pt.*y(1).^3.*y(2) - ...
4.*dK1t.*K2t.^3.*Pt.*y(1).^3 + 2.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^3 - 14.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2).^2 + 6.*dK1t.*K2t.^3.*Pt.*y(1).^2.*y(2) + ...
6.*dK1t.*K2t.^3.*Pt.*y(1).^2 + 8.*dK1t.*K2t.^3.*Pt.*y(1).*y(2).^2 - 16.*dK1t.*K2t.^3.*Pt.*y(1).*y(2) + 8.*dK1t.*K2t.^3.*Pt.*y(1) + 8.*dK1t.*K2t.^3.*Pt.*y(2) - ...
8.*dK1t.*K2t.^3.*Pt + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^3 + 16.*dK1t.*K2t.*Pt.^2.*y(1).^4.*y(2).^2 - 32.*dK1t.*K2t.*Pt.^2.*y(1).^3.*y(2).^3 + ...
16.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^3 - 48.*dK1t.*K2t.*Pt.^2.*y(1).^2.*y(2).^2 + 32.*dK1t.*K2t.*Pt.^2.*y(1).*y(2).^2))/(K1t.*K2t.*Pt.*y(1).*(- K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^4 + ...
2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^3.*y(2) + K1t.^2.*K2t.^2.*Pt.*y(1).^3 - 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^3 + ...
6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2).^2 - 6.*K1t.^2.*K2t.^2.*Pt.*y(1).^2.*y(2) + 2.*K1t.^2.*K2t.^2.*Pt.*y(1).^2 - 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^4 + ...
24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^3 - 24.*K1t.^2.*Pt.^2.*y(1).^3.*y(2).^2 + 8.*K1t.^2.*Pt.^2.*y(1).^3.*y(2) + 8.*K2t.^2.*y(1).^3.*y(2).^3 + ...
8.*K2t.^2.*y(1).^3.*y(2).^2 - 8.*K2t.^2.*y(1).^3.*y(2) - 8.*K2t.^2.*y(1).^3 - 8.*K2t.^2.*y(1).^2.*y(2).^3 + 24.*K2t.^2.*y(1).^2.*y(2).^2 + ...
8.*K2t.^2.*y(1).^2.*y(2) - 24.*K2t.^2.*y(1).^2 - 32.*K2t.^2.*y(1).*y(2).^2 + 32.*K2t.^2.*y(1).*y(2) - 32.*K2t.^2.*y(2) + 32.*K2t.^2 + 48.*Pt.*y(1).^3.*y(2).^3 + ...
32.*Pt.*y(1).^3.*y(2).^2 - 16.*Pt.*y(1).^3.*y(2) - 48.*Pt.*y(1).^2.*y(2).^3 + 32.*Pt.*y(1).^2.*y(2).^2 - 48.*Pt.*y(1).^2.*y(2) - 64.*Pt.*y(1).*y(2).^2 + 64.*Pt.*y(2)));
end
function Y = runge_kutta_RK4(f,T,Y0)
N = numel(T);
n = numel(Y0);
Y = zeros(N,n);
Y(1,:) = Y0;
for i = 2:N
t = T(i-1);
y = Y(i-1,:);
h = T(i) - T(i-1);
k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
end
end
Jason Louison
on 22 Dec 2022
Edited: Jason Louison
on 22 Dec 2022
As said, already the changes made to the RK integration scheme made your code wrong:
k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
is not the same as
k0 = h*f(t,y)';
k1 = h*f(t,y+k0*0.5)';
k2 = h*f(t,y+k1*0.5)';
k3 = h*f(t,y+k2)';
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
And the time-dependent parameters are incorporated in f:
dK1t = sin(3.14*t);
dK2t = 2*cos(3.14*t);
K1t = 60;
K2t = 25;
dPt = 1;
Pt = 1;
Jason Louison
on 22 Dec 2022
Don't you see that you update Y as
Y(i,:) = y + h^2/6*(k0+2*k1+2*k2+k3);
instead of
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
when the ki were defined in the usual way as
k0 = f(t,y);
k1 = f(t+0.5*h,y+k0*0.5*h);
k2 = f(t+0.5*h,y+k1*0.5*h);
k3 = f(t+h,y+k2*h);
?
For solving differential equations that are not dependent on time explicity, you would do the following:
Say we have a differential equation dy/dt = f(y) = 1 + y. Since the variable dy/dt itself does not depend on time, but on y, you would do:
k1 = f(y(t))
k1 is just the differential equation evaluated at a given value of y. t is the time, but y is NOT dependent on time explicitly. It is to express the value of y at a given time t, y(t).
k2 = f(y(t)+0.5*h*k1)
k2 is the estimation of dy/dt = f(y) at a half of a time step forward. from f(y). Since k1 is the derivative of y, you can essentially approximate y(t+0.5h) by doing this. k3 takes it a step further by using the esitmated derivative k2 and repeating the process:
k3 = f(y(n)+0.5*h*k2)
k4 is the final iteration in this process. It estimates the value of dy/dt at one full time step away from f(y) by usin k3 as the derivate.
k4 = f(y+h*k3)
to get the approximation of y at the next time step, you use this formula:
y(n+1) = y(n)+ (h/6)*(k1*2k2*2k3+k4)
After obtaining this new value of y, the process is repeated. This is the method I use for numerically integrating functions that are not explicity a function of time. Here is an example code:
%Sprung Mass Damper Simulation%
k = 1000; %Spring Constant
b = 600; %Damping Constant
m = 1000; %Mass
g = 9.81; %Gravity
p = k/m;
q = b/m;
tf = 120;
N = 10000;
h = tf/N;
t = linspace(0,tf,N+1);
x = zeros(1,N+1);
v = zeros(1,N+1);
x(1) = 0;
v(1) = 0;
for n=1:N %fourth-order runge kutta method algorithm
k1a = h*f(x(n),v(n),p,q,g);
k1v = h*v(n);
k2a = h*f(x(n)+1/2*k1v,v(n)+1/2*k1a,p,q,g);
k2v = h*(v(n)+1/2*k1a);
k3a = h*f(x(n)+1/2*k2v,v(n)+1/2*k2a,p,q,g);
k3v = h*(v(n)+1/2*k2a);
k4a = h*f(x(n)+k3v,v(n)+k3a,p,q,g);
k4v = h*(v(n)+k3a);
v(n+1) = v(n)+(1/6)*(k1a+2*(k2a+k3a)+k4a);
x(n+1) = x(n)+(1/6)*(k1v+2*(k2v+k3v)+k4v);
end
a = -p*x-q*v-g;
% It was hard to see the displacement when all three functions were graphed
% on the same axis, since x is so mauch smaller than v & a
% I split them onto separate lines
subplot(3,1,1);
plot(t,x)
title('displacement')
subplot(3,1,2)
plot(t,v)
title('velocity')
subplot(3,1,3)
plot(t,a) %plot displacement, velocity, and acceleration vs time
title('acceleration')
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function a = f(x,v,p,q,g)
a = -p*x-q*v-g;
end
Torsten
on 23 Dec 2022
Do you see the difference between the update of y in this code and the code this post is about ?
k1a = h*f(x(n),v(n),p,q,g);
k2a = h*f(x(n)+1/2*k1v,v(n)+1/2*k1a,p,q,g);
k3a = h*f(x(n)+1/2*k2v,v(n)+1/2*k2a,p,q,g);
k4a = h*f(x(n)+k3v,v(n)+k3a,p,q,g);
v(n+1) = v(n)+(1/6)*(k1a+2*(k2a+k3a)+k4a);
Here, you use 1/6 as multiplier.
k0 = h*f(t,y)';
k1 = h*f(t,y+k0*0.5)';
k2 = h*f(t,y+k1*0.5)';
k3 = h*f(t,y+k2)';
Y(i,:) = y + h/6*(k0+2*k1+2*k2+k3);
Here, you use h/6 as multiplier.
Jason Louison
on 23 Dec 2022
K1 = f(y(n)) K2 = f(y(n)+0.5*h*K1)) K3 = f(y(n)+0.5*h*K2)) K4 = f(y(n)+h*K3)
y(n+1) = (h/6)*(K1+2*(K2+K3)+K4)
instead of
K1 = f(y(n)) K2 = f(y(n)+0.5*K1)) K3 = f(y(n)+0.5*K2)) K4 = f(y(n)+K3)
y(n+1) = (h/6)*(K1+2*(K2+K3)+K4)
It was a mistype in the chemical equilibrium code, not in the spring code.
And correcting it will lead to unphysical results.
Jason Louison
on 23 Dec 2022
Jason Louison
on 24 Dec 2022
Edited: Jason Louison
on 24 Dec 2022
Torsten,
You are right, this is very strange. The strangest part is that I get normal, correct results when I write the algorithm as h^2/6*(k1+2*k2+2*k3+k4), but if I try writing the algorithm properly it, it breaks the simulation. It doesn't make any sense, as my other simulations work correctly with the algorithm I just explained (With the picture of the algorithm). Do you know what might be causing this? Attached below is the results of the "broken" simulation.
I will try re-setting the code back to the way you originally did it, but I am experienced in doing these simulations on spreadsheets, not code, and that is why I modified it the way I did. On a spreadsheet, using the runge kutta algorithm, you cannot include time in the alogorithm, even if a function changes with time. It would already factor it in. I did it this way in my spring damper system code as well, and it worked. Maybe because parameters are changing with time, I have to include time in the function. Anyways, I will get back to you soon with results of the new code.
More Answers (0)
Categories
Find more on Chemistry in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
