How to calculate dx/dt in Matlab

I have the following equation , where DOT is a vector [50;51....99], DOT^* is a constant = 100, Time has the same length as DOT

I know th symbolic solution which is for that:

But I need to get Kla value numerically, how do I do that in Matlab, I tried

kla = (gradient(DOT) ) ./ (gradient(Time) .* (DOT_star - DOT));

but I get very different answer from the analytical one. Any Idea what can be wrong in my equation ?

[Torsten]

You say DOT = [50;51;...;99].

And what is the corresponding vector "Time" ?

[Walter Roberson]

I am not sure where your DOT_(end) comes from?

syms DOT__star DOT(t) K_La DOT_t1 t1
dDOT = diff(DOT,t)

dDOT(t) = 

eqn = dDOT == (DOT__star - DOT)*K_La

eqn(t) = 

ic = DOT(t1) == DOT_t1

ic = 

sol = dsolve(eqn, ic)

sol = 

solve(sol, K_La)

ans = 

Also instead of gradient(DOT)./gradient(time) I suggest you gradient(DOT, time)

[Torsten]

solve(sol == DOT, K_La)

instead of

solve(sol, K_La)

Since the resulting equation for K_la holds for every t, it also holds for t = t(end).

(All under the assumption that K_la is constant).

[Walter Roberson]

Ah,

syms DOT__star DOT(t) K_La DOT_t1 t1
dDOT = diff(DOT,t)

dDOT(t) = 

eqn = dDOT == (DOT__star - DOT)*K_La

eqn(t) = 

ic = DOT(t1) == DOT_t1

ic = 

sol = dsolve(eqn, ic)

sol = 

solve(sol == DOT, K_La)

ans = 

[Sami]

thanks very much for the analytical solution, any idea how to do the numerical one to get Kla value ?

[Torsten]

You did it correctly, but - as I noted - if

kla = (gradient(DOT) )  ./ (gradient(Time) .* (DOT_star - DOT));

does not return an almost constant vector, your data vectors (Time,DOT) do not follow the mass transfer low

d(DOT)/dt = kla*(DOT_star - DOT)

If you want the "best" kla - value for your data, use an optimizer (e.g. lsqnonlin).

Here is some code:

Time = ...;
DOT = ...;
DOT_star = ...;
fun = @(kla) DOT - DOT_star - (DOT(1) - DOT_star).*exp(-kla*(Time-Time(1)));
kla = lsqnonlin(fun,1)

[Sami]

Thanks, when can we consider dDOT/dt equals gradient(DOT) ./ gradient(Time) ?

[Torsten]

Approximately: always. But I don't understand the sense of your question.

[Sami]

since you said :

------------------

You did it correctly, but - as I noted - if ......

does not return an almost constant vector, your data vectors (Time,DOT) do not follow the mass transfer low

------------------

so I was wondering how can one check when dDOT/dt can be considered equal to gradient(DOT) ./ gradient(Time) , and when not ?

[Torsten]

The easiest way is if you include your Time/DOT data as a two-column matrix.

Or did you just invent such vectors and now wonder why they don't show the behaviour of the model equation ?

Because your DOT vector above ([50;51....99]) surely does not follow the model equation (at least for an equidistant corresponding vector Time).

[Sami]

I do have two column matrix (DOT, TIME) with equadistant time, still the answer is wong,

I put here some values:

Time = [

15.5189990000000 15.5194990000000 15.5199990000000 15.5204990000000 15.5209990000000 15.5214990000000 15.5219990000000 15.5224990000000 15.5229990000000 15.5234990000000 15.5239990000000 15.5244990000000 15.5249990000000 15.5254990000000 15.5259990000000 15.5264990000000 15.5269990000000 15.5274990000000 15.5279990000000 15.5284990000000 15.5289990000000 15.5294990000000 15.5299990000000 15.5304990000000 15.5309990000000 15.5314990000000 15.5319990000000 15.5324990000000 15.5329990000000 15.5334990000000 15.5339990000000 15.5344990000000 15.5349990000000 15.5354990000000 15.5359990000000 15.5364990000000 15.5369990000000 15.5374990000000 15.5379990000000 15.5384990000000 15.5389990000000 15.5394990000000 15.5399990000000 15.5404990000000 15.5409990000000 15.5414990000000 15.5419990000000]

DOT = [86.2934346288157 88.3330265127609 88.5384902983721 86.2715655041141 86.2499514074551 87.9439995156669 93.3559053355719 100.588355429879 105.134156185296 105.998293130410 101.914607317850 97.2442257193124 96.1705619186325 97.6457751783728 102.639322947358 106.338201652074 103.282744164835 99.3874791803122 98.3360407005509 98.1827495464303 98.4462091917985 98.4677691127529 98.6027332770350 99.3776569743805 101.164673911659 102.507302979436 101.756315732742 100.862114099507 100.796601964307 100.746832837656 100.715064263002 100.703558708876 100.714578643816 100.750386536354 100.813244855022 100.905416068358 101.029162644892 101.186747053159 101.341304721314 101.205048020864 100.712709159760 100.334695799161 100.192106610432 100.121980159814 100.129391137103 100.219415163613 100.329039499815]

Dot_star = 100;

I get different Kla for the analytical and the numerical solutions, I am kind of double checked everything, still not seeing where the error might come from!!

[Torsten]

This is the best approximation to your data with a resulting kla value of 389.

I think you will admit that with such swinging data, taking numerical gradients leads to nowhere.

Time = [15.5189990000000	15.5194990000000	15.5199990000000	15.5204990000000	15.5209990000000	15.5214990000000	15.5219990000000	15.5224990000000	15.5229990000000	15.5234990000000	15.5239990000000	15.5244990000000	15.5249990000000	15.5254990000000	15.5259990000000	15.5264990000000	15.5269990000000	15.5274990000000	15.5279990000000	15.5284990000000	15.5289990000000	15.5294990000000	15.5299990000000	15.5304990000000	15.5309990000000	15.5314990000000	15.5319990000000	15.5324990000000	15.5329990000000	15.5334990000000	15.5339990000000	15.5344990000000	15.5349990000000	15.5354990000000	15.5359990000000	15.5364990000000	15.5369990000000	15.5374990000000	15.5379990000000	15.5384990000000	15.5389990000000	15.5394990000000	15.5399990000000	15.5404990000000	15.5409990000000	15.5414990000000	15.5419990000000];
DOT = [86.2934346288157	88.3330265127609	88.5384902983721	86.2715655041141	86.2499514074551	87.9439995156669	93.3559053355719	100.588355429879	105.134156185296	105.998293130410	101.914607317850	97.2442257193124	96.1705619186325	97.6457751783728	102.639322947358	106.338201652074	103.282744164835	99.3874791803122	98.3360407005509	98.1827495464303	98.4462091917985	98.4677691127529	98.6027332770350	99.3776569743805	101.164673911659	102.507302979436	101.756315732742	100.862114099507	100.796601964307	100.746832837656	100.715064263002	100.703558708876	100.714578643816	100.750386536354	100.813244855022	100.905416068358	101.029162644892	101.186747053159	101.341304721314	101.205048020864	100.712709159760	100.334695799161	100.192106610432	100.121980159814	100.129391137103	100.219415163613	100.329039499815];
DOT_star = 100;
fun = @(kla) DOT - DOT_star - (DOT(1) - DOT_star).*exp(-kla*(Time-Time(1)));
kla = lsqnonlin(fun,1)

Local minimum possible.

lsqnonlin stopped because the final change in the sum of squares relative to 
its initial value is less than the value of the function tolerance.

kla = 388.7042

plot(Time,DOT)
hold on
plot(Time,DOT_star + (DOT(1) - DOT_star).*exp(-kla*(Time-Time(1))))

Answer 1

[Torsten]

If your results from the theoretical equation for K_la and from the numerical formula don't correspond, it's indicating that the mass transfer law

du/dt = K_la * (u* - u)

doesn't apply for your data.

Note that usually the formula

kla = (gradient(DOT) )  ./ (gradient(Time) .* (DOT_star - DOT));

shouldn't work at all because gradient(DOT) and DOT don't have the same number of elements.

[Sami]

actually, in this case both does, gradient gives the ame vector length of DOT