ODE45 and dsolve result difference

Hi there! I am tring to solve a system of differential equations using both ode45 and dsolve. However, the outputs graph of both methods are very different.

Here's the ode45's code:

clear

clc

% Parameter

k1 = 1000000;

k2 = 10;

k3 = 1000000;

k4 = 10;

A = 0.00001;

B = 0.00002;

alpha = 1.5;

% Initial Condition and Time Range

y0 = [0; 0; 0; 1 ]; % Initial Condition

tspan = [0 0.5]; % Time Range

% Solving System of DE

[t, y] = ode45(@(t,y) [k1*A*alpha*y(3) - k2*y(1) + k3*alpha*B*y(2) - k4*y(1);

k1*A*y(4) - k2*y(2) - k3*alpha*B*y(2) + k4*y(1);

-k1*A*alpha*y(3) + k2*y(1) + k3*B*y(4) - k4*y(3);

-k1*A*y(4) + k2*y(2) - k3*B*y(4) + k4*y(3)], tspan, y0);

Then the dsolve code:

syms v(t) y(t) w(t) z(t)

% Parameter

k1 = 1000000;

k2 = 10;

k3 = 1000000;

k4 = 10;

A = 0.00001;

B = 0.00002;

alpha = 1.5;

% Defining System of Differential Equations

eqns = [

diff(y) == k1*A*alpha*z - k2*y + k3*alpha*B*v - k4*y,

diff(v) == k1*A*w - k2*v - k3*alpha*B*v + k4*y,

diff(w) == -k1*A*alpha*z + k2*y + k3*B*w - k4*z,

diff(z) == -k1*A*w + k2*v - k3*B*w + k4*z

];

initCond = [y(0) == 0, v(0) == 0, w(0) == 1, z(0) == 0];

% Solving DEs

[vSol(t), ySol(t), wSol(t), zSol(t)] = dsolve(eqns,initCond);

% Showing Results

v = vpa(vSol(t))

y = vpa(ySol(t))

w = vpa(wSol(t))

z = vpa(zSol(t))

% Plotting the graph

fplot(v, [0, 0.5], 'LineWidth', 2);

hold on;

fplot(y, [0, 0.5], 'LineWidth', 2);

fplot(w, [0, 0.5], 'LineWidth', 2);

fplot(z, [0, 0.5], 'LineWidth', 2);

xlabel('Time (t)');

ylabel('Value');

legend('v(t)', 'y(t)', 'w(t)', 'z(t)');

hold off;

For clarification, y(1), y(2), y(3) and y(4) correspond to, respetively, y, v, z, w.

Graph-wise. ode45 and dsolve's graph are, respectively,

and .

Answer 1

[Sam Chak]

The order of the state variables between and is swapped in the symbolic approach.

syms v(t) y(t) w(t) z(t)

% Parameter
k1 = 1000000;
k2 = 10;
k3 = 1000000;
k4 = 10;
A = 0.00001;
B = 0.00002;
alpha = 1.5;

% Defining System of Differential Equations
eqns = [
    diff(y) ==  k1*A*alpha*w - k2*y + k3*alpha*B*v - k4*y,
    diff(v) ==  k1*A*z - k2*v - k3*alpha*B*v + k4*y,
    diff(w) == -k1*A*alpha*w + k2*y + k3*B*z - k4*w,
    diff(z) == -k1*A*z + k2*v - k3*B*z + k4*w
];

%  [diff(y) ==  k1*A*alpha*y(3) - k2*y(1) + k3*alpha*B*y(2) - k4*y(1);
%   diff(v) ==  k1*A*y(4) - k2*y(2) - k3*alpha*B*y(2) + k4*y(1);
%   diff(w) == -k1*A*alpha*y(3) + k2*y(1) + k3*B*y(4) - k4*y(3);
%   diff(z) == -k1*A*y(4) + k2*y(2) - k3*B*y(4) + k4*y(3)];

initCond = [y(0) == 0, v(0) == 0, w(0) == 0, z(0) == 1];

% Solving DEs
[vSol(t), ySol(t), wSol(t), zSol(t)] = dsolve(eqns,initCond);

% Showing Results
v = vpa(vSol(t));
y = vpa(ySol(t));
w = vpa(wSol(t));
z = vpa(zSol(t));

% Plotting the graph
hold on;
fplot(w, [0, 0.5], 'LineWidth', 2);
fplot(v, [0, 0.5], 'LineWidth', 2);
fplot(y, [0, 0.5], 'LineWidth', 2); 
fplot(z, [0, 0.5], 'LineWidth', 2);
xlabel('Time (t)');
ylabel('Value');
legend('w(t)', 'v(t)', 'y(t)', 'z(t)');
hold off; grid on