Solving a non-linear ODE with coupled algebraic equations
5 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Ali Muhammad Hassaan
am 29 Aug. 2020
Kommentiert: Alan Stevens
am 30 Aug. 2020
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?
17 Kommentare
Akzeptierte Antwort
Alan Stevens
am 29 Aug. 2020
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
8 Kommentare
Alan Stevens
am 30 Aug. 2020
"...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
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!