Code Won't Finish running

3 次查看（过去 30 天）
Will Jeter
Will Jeter 2020-11-3
回答： Mathieu NOE 2020-11-3
Cant get my ODEs to run correctly with my X(2) value. Code will run but it doesn't complete and get a giant debugging page if I pause. Please help
[t,X] = ode45(@F,[0 340],[0; 400]);
figure
plot(t,X(:,1),'o')
legend('X')
xlabel('Reaction time (min)')
ylabel('X')
% Find t@X = 0.9
indx = find(abs(X(:,1)-0.9) < 0.01)
format
fprintf('Time (min) to achieve X = 0.9 is %f', t(indx))
fprintf('Conversion is %f', X(indx,1))
function dCdt = F(t,X)
R = 8.314; % J/mol/K
CA0 = 0.5; % M
CB0 = 0.75; % M
Cp = 3.8; % J/g/K
p = 0.75; %g/cm^3
dHrxn = -145000; % J/mol
V = 50; % L
NA0 = V*CA0; % mol
T0 = 400; %K
UA = 100; % W/K
ka = 1.4*10^7*exp(-7700/X(2)); % L/mol/min
dCdt = zeros(2,1);
dCdt(1,1) = ka*CA0*(1-X(1))*(1.5-X(1));
dCdt(2,1) = (-dHrxn/(p*Cp))*CA0^2*ka*(((1-X(1))*(1.5-X(1)))/1+(X(1)));
end

请先登录，再回答此问题。

回答（2 个）

Alan Stevens
Alan Stevens 2020-11-3
Use ode15s instead of ode45. Also, you might as well set the timespan to be [0 0.1] as nothing much changes after that!
Mathieu NOE
Mathieu NOE 2020-11-3
hello
your code works. I first reduced your time span because it was taking forever...
I also change the ode solver because - according to the plot - you shoud better use a solver for stiff ode.
speed up the whole thing quite significantly
the time transition you are looking for is around t = 0.02 s so no need to create a huge time span up to 340 s - unless there is something special happening at the end of the simulation ?
code more or less unchanged - minor stuff :
% [t,X] = ode45(@F,[0 340],[0; 400]);
[t,X] = ode15s(@F,[0 0.1],[0; 400]);
figure(1),
plot(t,X(:,1),'o')
legend('X')
xlabel('Reaction time (min)')
ylabel('X')
% Find t@X = 0.9
indx = find(abs(X(:,1)-0.9) < 0.01)
format
% fprintf('Time (min) to achieve X = 0.9 is %f', t(indx))
% fprintf('Conversion is %f', X(indx,1))
disp(['Time (min) to achieve X = 0.9 is : ' num2str(t(indx)) ' s']);
disp(['Conversion is : ' num2str(X(indx,1)) ]);
function dCdt = F(t,X)
R = 8.314; % J/mol/K
CA0 = 0.5; % M
CB0 = 0.75; % M
Cp = 3.8; % J/g/K
p = 0.75; %g/cm^3
dHrxn = -145000; % J/mol
V = 50; % L
NA0 = V*CA0; % mol
T0 = 400; %K
UA = 100; % W/K
ka = 1.4*10^7*exp(-7700/X(2)); % L/mol/min
dCdt = zeros(2,1);
dCdt(1,1) = ka*CA0*(1-X(1))*(1.5-X(1));
dCdt(2,1) = (-dHrxn/(p*Cp))*CA0^2*ka*(((1-X(1))*(1.5-X(1)))/1+(X(1)));
% dCdt = [ka*CA0*(1-X(1))*(1.5-X(1));(-dHrxn/(p*Cp))*CA0^2*ka*(((1-X(1))*(1.5-X(1)))/1+(X(1)));];
% % do the same as original code (the 3 line above are equivalent to this single line)
end

类别

在 Help Center 和 File Exchange 中查找有关 Watermarking 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by