Info

此问题已关闭。 请重新打开它进行编辑或回答。

Why does defining a system of ODEs a different, but similar way yield very different solutions?

1 次查看（过去 30 天）
Riley Stein
Riley Stein 2020-8-13
关闭： MATLAB Answer Bot 2021-8-20
I have a system of ODEs that I am trying solve, however the correct solutions are not being returned (I have another script to which I'm comparing the produced values). As I am trying to figure out why, I defined the function of ODEs as a column vector, instead of how it was originally defined as row vectors, and got a different answer.
My original code produces values closest in magnitude to what I want and is defined as follows:
Enfcn = @(t) Ei + v * t;
kred = @(t) 225000 * exp(-8) * exp(-(a * f) * (Enfcn(t) - Eo));
kox = @(t) 225000 * exp(-8) * exp(((1-a) * f) * (Enfcn(t) - Eo));
C0 = [0.1, 0, 0, 0, 0, 1e6, 0];
tspan = 0:0.01:15;
[Time, Concentrations] = ode15s(@(t, C) mechanism(t, C, kred, kox), tspan, C0);
function dC = mechanism(t, C, kred, kox)
%Rate Constants (ish)
k1 = 130;
k2 = 1;
pKa = 2.5;
%Protonation Rates
ka = 1e1;
Ka = 10^(-pKa) * 1e9;
kb = Ka * ka;
% Variables
CNiII = C(1);
CNiIISH = C(2);
CNiISH = C(3);
CNiIIIH = C(4);
CNiIIH = C(5);
Ce = C(6); %I know this is unused, but it doesn't hurt having it here
CH = C(7);
CH2 = C(8); %Same with this
% Rate Laws / ODEs
dNiII = -ka*CNiII*CH + kb*CNiIISH + k2*CNiIIH*CH;
dNiIISH = -kb*CNiIISH + ka*CNiII*CH - kred(t)*CNiIISH + kox(t)*CNiISH;
dNiISH = kred(t)*CNiIISH - kox(t)*CNiISH - k1*CNiISH;
dNiIIIH = -kred(t)*(CNiIISH + CNiIIIH) + kox(t)*(CNiISH + CNiIIH);
dNiIIH = k1*CNiISH - kred(t)*CNiIIIH + kox(t)*CNiIIH;
de = kred(t)*CNiIIIH - kox(t)*CNiIIH - k2*CNiIIH*CH;
dH = -k2*CNiIIH*CH - ka*CNiII*CH + kb*CNiIISH;
dH2 = k2*CNiIIH*CH;
% Assign Output Variables
dC(1,:) = dNiII;
dC(2,:) = dNiIISH;
dC(3,:) = dNiISH;
dC(4,:) = dNiIIIH;
dC(5,:) = dNiIIH;
dC(6,:) = de*(96485.33289/1e9); %This needs to multiplied at some point by this conversion, I do it when I plot, I just added it here
dC(7,:) = dH;
dC(8,:) = dH2;
When changing it to this, I gain the plot shape that I want, but not the correct values (in fact this one is several orders of magnitude off), and the integration tolerance of ode15s can't be met even if I set RelTol and AbsTol to values from 1e-12~1e-20:
dC = [-ka*CNiII*CH + kb*CNiIISH + k2*CNiIIH*CH;
-kb*CNiIISH + ka*CNiII*CH - kred(t)*CNiIISH + kox(t)*CNiISH;
kred(t)*CNiIISH - kox(t)*CNiISH - k1*CNiISH;
-kred(t)*(CNiIISH + CNiIIIH) + kox(t)*(CNiISH + CNiIIH);
k1*CNiISH - kred(t)*CNiIIIH + kox(t)*CNiIIH;
kred(t)*CNiIIIH - kox(t)*CNiIIH - k2*CNiIIH*CH;
-k2*CNiIIH*CH - ka*CNiII*CH + kb*CNiIISH;
k2*CNiIIH*CH];
I am confused as to why these give different results when from my understanding they were synonymous, but I am also still new to this. Any insight would be very appreciated. Thanks in advance.

此问题已关闭。

回答（1 个）

James Tursa
James Tursa 2020-8-13
Do you simply need to apply that factor in the 6th spot?
(kred(t)*CNiIIIH - kox(t)*CNiIIH - k2*CNiIIH*CH)*(96485.33289/1e9);
  1 个评论
Riley Stein
Riley Stein 2020-8-13
编辑：Riley Stein 2020-8-13
I figured out the problem because of your comment. The original function can be in any order because I later define what is exaclty what, whereas the column vector is sensitive to exactly where the ODE is placed, which makes sense. I'm sorry for such a senseless error. I was too excited about getting the plot shape correct that I didn't even consider the actual variable naming. Although I am still perplexed by why the ODEs are not being solved the same.

此问题已关闭。

标签

Community Treasure Hunt

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

Start Hunting!

Translated by