How can I solve the fsolve function inside an ode15s function?

1 view (last 30 days)
Hi,
I need to solve an fsolve function inside an ode15s function. However, the dependent variable in the differential equation is contained in the system of linear equations. How can I solve this? This is what I have but it says there is a problem:
%file1.m
clc; clear all; close all; format compact;
global A_T Fao P_tot DELTAHr E A R To Tf
%Datos
A_T = pi*1^2;%m
Fao = (66000*1000)/92.14; %kmol/año
P_tot = 50; %bar
DELTAHr = 49974; %KJ/Kmol
%Cálculo de k
E = 148114;
A = 5.73E6;
R = 8.3144;
To = 600+273;
Tf = 635+273;
Zo = 0;
Zf = 100;
yo = [To, 0];
[F fval flag] = fsolve('fun',[3 10])
[T, xa] = ode15s('odefun', [Zo Zf], yo)
%file2.m
function F = fun(var)
global A_T Fao P_tot DELTAHr E A R To Tf FA FR CPA CPB CPR CPS CQP
FB = var(1)
FS = var(2)
Q = var(3)
F(1) =-FB + Fao*(5-x_a) + 0.9*Q ;
F(2) =-FS + Fao*(5/9 + x_a) + 0.1*Q;
F(3) = -Q +(FA*CPA*(T-To)+ FB*CPB*(T-To)+ FR*CPR*(T-To)+ FS*CPS*(T-To))/(CQP*(To-298));
end
%file3.m
function dvdz = odefun(T, x_a)
global A_T Fao P_tot DELTAHr E A R To CPA CPB CPR CPS CQP FA FR
K = A*exp(-E/R*T);
CPA = 0.29 + 47.052E3*T - 15.716E6*T^2;
CPB = 3.249 + 0.422E3*T + 0.083E-5*T^-2;
CPR = -.206 + 39.064E3*T - 13.301E6*T^2;
CPS = 1.702 + 9.081E3*T - 2.164E6*T^2;
CQP = CPB*0.9 + CPS*0.1;
FA = Fao*(1-x_a);
FR = Fao*x_a;
var0 = [100 100 100];
[var] = fsolve('fun', var0);
var(1) = FB ;
var(2) = FS ;
var(3) = Q;
dvdz(1) = (A_T/Fao)*K*(P_tot/R*T)^1.5*(Fao*(1-x_a))*(((Fao*(5-x_a) + 0.9*Q))^0.5/(6.55*Fao + Q)^1.5);
dvdz(2) = (Fao*DELTAHr/(FA*CPA*(T-To)+ FB*CPB*(T-To)+ FR*CPR*(T-To)+ FS*CPS*(T-To)))*dvdz(1);
dvdz = dvdz'
end
Thank you!
  1 Comment
Torsten
Torsten on 18 May 2015
From your code above, I don't understand which system you are trying to solve.
If it helps: Usually, there is no Need to solve algebraic equations separately. The ODE solvers are suited to solve a mixture of algebraic and differential equations. Take a look at the differential-algebraic example under
Best wishes
Torsten.

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!