Solving a non-linear ODE with coupled algebraic equations
10 次查看(过去 30 天)
显示 更早的评论
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 个评论
采纳的回答
Alan Stevens
2020-8-29
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 个评论
Alan Stevens
2020-8-30
"...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
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!