Solving a non-linear ODE with coupled algebraic equations

I'm trying to solve a differential equation which has two coupled algebraic equations and two independent variables.

, is the differential equation.

is the phase and varies from 0 to 1. The initial condition for the phase is

n and t* are defined as,

My knowledge of solving differential equations in matlab is limited and I'm unable to couple the equations with the differential equation.

Also, what would be the best way to solve this problem, should I use the standard ode45 or use DAE solver?

[Alan Stevens]

Is p2 being differentiated with respect to t or t*?

What is the initial Temperature?

What is the timespan of the integration?

[Ali Muhammad Hassaan]

Sorry, the initial tempperature is T(0)=0 and the timespan is 400 s.

[Alan Stevens]

What about t and t* ?

If you ae integrating wrt t not t*, then you need something else to relate t and T.

If wrt t* then T(0) is not consistent with t* = 0.

[Ali Muhammad Hassaan]

The temperature T is a function of time and has the same timespan of 100s. The initial temperature is T(0)=0. Maximum temperature T varies from 0 to 500 C.

[Alan Stevens]

So your first equation means dp/dt* = ...?

If so, what happens at t* = 0? It implies dp/dt* is infinite or undefined (divide by zero). Also, The relation between T and t* is inconsistent.

[Ali Muhammad Hassaan]

When T is 0, t* is inf, and dp/dt becomes 0.

[Alan Stevens]

So, you aren't integrating with respect to t* then?

Also, if p2(0) = 0, and dp2dt(0) is zero then there is nothing that allows p2 to change - it will stay at zero!

[Ali Muhammad Hassaan]

Fair enough, say, if I start at room temperature then the time evolution will not be equal to zero and the temperature profile as well as dp/dt evolve with time?

Sorry for the confusion.

[Alan Stevens]

If you start at room temperature then, according to your relation between t* and T, the initial value of t* will be around 600.

And if t* is finite, then with p2(0) = 0, dp2/dt will also be zero, so p2 will still not change!

[Ali Muhammad Hassaan]

Yeah, initially, the phase p2 does not evolve with time and does so only after about 20 seconds. Actually, I am trying to replicate results from a paper and either the information cited is incorrect or there is something that I am missing.

[Alan Stevens]

There is still a problem at 20 secs! If p2 is zero, dp2/dt is still either zero or infinite (depending on p being on numerator or denominator - the latter determined by the value of n), neither of which is very useful for integrating!

Do you have a copy of the paper?

[Ali Muhammad Hassaan]

Yeah, how do I send it to you?

[Alan Stevens]

Use the attachments symbol (the paper-clip symbol) to upload it here. If it's not in electronic form, either scan it first, or find a link to an (open) on-line version.

There's no guarantee I'll understand it, but I'll have a look!

[Ali Muhammad Hassaan]

Here's the paper. Take a look at section 4.2.1.2

[Alan Stevens]

Ok. Leave it with me. It might take me some time to read and understand - if I understand it at all!

[Ali Muhammad Hassaan]

Okay, I will wait for your answer. Thanks for taking the time to go through it.

[Ali Muhammad Hassaan]

Hey Alan, is it possible to use the date from the graph for the temperature against time and use that data for the values of T (temperature)?

Answer 1

[Alan Stevens]

Yes, that's exactly what I've done in the program below. However, in order to get anything like the results in the paper (which I don't understand!) I've had to multiply the expression for dp2/dt by a factor of 0.01 (there might be an issue with units that require this?). Even so, the results (see below) don'tlook exactly like those of the paper, but it's the best I can do.

% Basic data
a = 0.394;
b = 0.107;
Qs = 30;  %kJ/mol
Qd = 130; %kJ/mol
tr1 = 600; %s
Tr1 = 375 + 273.15; %K
%T = 300 + 273.15; %K
R = 8.3145; %j/(mol.K)

% Timespan and initial condition
tspan = [0 400];
p20 = 10^-8;

% Call ode function
[t, p2] = ode15s(@fn,tspan,p20,[],a,b);

% Calculate temperatures as a function of time
T = min(300/20*t, 300) + 273.15;

% Plot results 
figure
plot(t, p2,'--',t,1-p2,'b:'), grid
xlabel('Time (s)'),ylabel('p2')
yyaxis right
plot(t,T-273.15,'k')
axis([0 400 0 500])
ylabel('Temperature (C)')
legend('p2', 'p1', 'Temperature')

% dp2/dt function
function dp2dt = fn(t, p2, a,b)

        tstar = tstarfn(t);
        n = 0.5 - a*p2^b;
        if t<20
            dp2dt = 0;
        else
            dp2dt = 0.01*n*p2^(1-1/n)/tstar; % Note: Multiplied by extra factor of 0.01 to get better agreement with paper
        end


end

function tstar = tstarfn(t)
        
        Qs = 30;  %kJ/mol
        Qd = 130; %kJ/mol
        tr1 = 600;
        Tr1 = 375 + 273.15; %K
        R = 8.3145; %j/(mol.K)

        T = min(300/20*t, 300) + 273.15;
        tstar= tr1*exp((Qs/(0.5*R) + Qd/R)*(1./T - 1/Tr1));
        
end

[Ali Muhammad Hassaan]

How about this, I extracted the data from the graph from the paper, here it is:

X = [0 ,0.5681,1.7043,2.2724,2.8405,3.4086,4.5448,5.1129,5.681, 6.8172,8.5215,9.0898,9.6577,10.7933,12.4982,14.2025,15.9068,18.1792,20.4516,22.1559,22.724,25,27.2688,29.5412,30.1093,31.8136,32.9498,34.086,50,75,100]

Y = [12.5, 25, 50, 62.5, 75, 100, 118.75, 125, 150, 175, 200, 212.5, 225, 250, 275, 300, 325, 350, 375, 400, 406.25, 425, 450, 468.75, 475, 487.5, 496.875, 500, 500, 500, 500]

This is the data from the second graph.

[Ali Muhammad Hassaan]

Actually, this is already more than I could do. I think I should accept the answer, please update if you have the time.

Thank you for your effort and time. Much appreciated!

[Ali Muhammad Hassaan]

Why use 0.5 instead of n in tstar, also the values of Qs and Qd are in KJ/mol, while R is in J/(mol.K). So if you put Qs = 30000 and Qd = 130000, you don't need the scaling factor anymore.

[Alan Stevens]

That's great! Well done. And you are quite right - it should be n not just 0.5 in tstar.

[Alan Stevens]

Even with your units changes there is still a mismatch. I've got better agreement now, with no need for an artificial delay of 20 seconds. However, I have had to reduce tr1 from 600 to 20 seconds.

% Basic data
a = 0.394;
b = 0.107;
Qs = 30;  %kJ/mol
Qd = 130; %kJ/mol
tr1 = 600; %s  
Tr1 = 375 + 273.15; %K
R = 8.3145; %j/(mol.K)

% Timespan and initial condition
tspan = [0 400];
p20 = 10^-8;      %%%%%%%%%%  Very small initial amount

% Call ode function
[t, p2] = ode15s(@fn,tspan,p20,[],a,b);

% Calculate temperatures as a function of time
T = min(300/20*t, 300) + 273.15;

% Plot results 
figure
plot(t, p2,'--',t,1-p2,'b:'), grid
xlabel('Time (s)'),ylabel('p2')
yyaxis right
plot(t,T-273.15,'k')
axis([0 400 0 500])
ylabel('Temperature (C)')
legend('p2', 'p1', 'Temperature')

% dp2/dt function
function dp2dt = fn(t, p2, a,b)
        
        n = 0.5 - a*p2^b;
        tstar = tstarfn(t,n);
        dp2dt = n*p2^(1-1/n)/tstar; 

end

function tstar = tstarfn(t,n)
        
        Qs = 30000;  %J/mol
        Qd = 130000; %J/mol
        tr1 = 600/30;  %s     %%%%%%%%%%%%%%%%%%%%% Modified to match paper results
        Tr1 = 375 + 273.15; %K
        R = 8.3145; %j/(mol.K)
       T = min(300/20*t, 300) + 273.15;
        tstar= tr1*exp( (Qs/(n*R) + Qd/R)*(1./T - 1/Tr1) );
        
end

I think I'll call it a day there though!

[Alan Stevens]

I get a pretty good match with the 500C case by using your temperature vs time relationship; again setting tr1 to be 20, but also modifying Tr1 to be 473C:

[Ali Muhammad Hassaan]

Hey Alan, going back to my initial question, how do you decide the type of solver to use? Also, can you share the code for the second graph so that I can compare the results and my approach. Again, thanks for your help.

[Alan Stevens]

"...how do you decide the type of solver to use? "

This is basically just an ODE, with some extra relationships (tstar vs T), so I tried ode45 first. In this case there were a couple of steps where ode45 produced some complex values, so I switched to ode15s, which turned out to be better.

Here is the code for the second problem:

% Basic data
a = 0.394;
b = 0.107;
Qs = 30;  %kJ/mol
Qd = 130; %kJ/mol
tr1 = 20; %s  
Tr1 = 375 + 273.15; %K
R = 8.3145; %j/(mol.K)

% Timespan and initial condition
tspan = [0 100];
p20 = 10^-8;      %%%%%%%%%%  Very small initial amount

% Call ode function
[t, p2] = ode15s(@fn,tspan,p20,[],a,b);

% Calculate temperatures as a function of time
T = Tmpfn(t);

% Plot results 
figure
plot(t, p2,'--',t,1-p2,'b:'), grid
xlabel('Time (s)'),ylabel('p2')
axis([0 100 0 1])
yyaxis right
plot(t,T-273.15,'k')
axis([0 100 0 500])
ylabel('Temperature (C)')
legend('p2', 'p1', 'Temperature')

% dp2/dt function
function dp2dt = fn(t, p2, a,b)
        
        n = 0.5 - a*p2^b;
        tstar = tstarfn(t,n);
        dp2dt = n*p2^(1-1/n)/tstar; 

end

function tstar = tstarfn(t,n)
        
        Qs = 30000;  %J/mol
        Qd = 130000; %J/mol
        tr1 = 20;  %s     %%%%%%%%%%%%%%%%%%%%% Modified to match paper results
        Tr1 = 473 + 273.15; %K  % Original was 375C  %%%%%%%%%%%%%%%%%%
        R = 8.3145; %j/(mol.K)
        T = Tmpfn(t);
        tstar= tr1*exp( (Qs/(n*R) + Qd/R)*(1./T - 1/Tr1) );
        
end

function T = Tmpfn(t)
     
     tm = [0 ,0.5681,1.7043,2.2724,2.8405,3.4086,4.5448,5.1129,5.681, 6.8172, ...
         8.5215,9.0898,9.6577,10.7933,12.4982,14.2025,15.9068,18.1792,20.4516,...
         22.1559,22.724,25,27.2688,29.5412,30.1093,31.8136,32.9498,34.086,50,75,100];
    Tmp = [12.5, 25, 50, 62.5, 75, 100, 118.75, 125, 150, 175, 200, 212.5, 225,...
        250, 275, 300, 325, 350, 375, 400, 406.25, 425, 450, 468.75, 475,...
         487.5, 496.875, 500, 500, 500, 500];
       T = interp1(tm, Tmp, t) + 273.15;
    
end