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