ODE45 and dsolve result difference

19 次查看（过去 30 天）
Bohan Li
Bohan Li 2024-3-30
编辑： Sam Chak 2024-3-30
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 .

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

采纳的回答

Sam Chak
Sam Chak 2024-3-30
编辑：Sam Chak 2024-3-30
Ran in:
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

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Symbolic Math Toolbox 的更多信息

标签

产品


版本

R2024a

Community Treasure Hunt

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

Start Hunting!

Translated by