MATLAB: Solving a non-linear ODE with coupled algebraic equations

non-linear odeode

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?

Best Answer

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

Best Answer

Related Question