conversion a symbolic function

2 次查看（过去 30 天）
bianchi fiammetta
bianchi fiammetta 2018-11-30
评论： Walter Roberson 2018-11-30
Good morning!
I do not manage to solve this error: I have to introduce a symbolic function in another function (not symbolic). I tried to convert by vpa and double, but the program does not work.
Thanks you
clear all
clc
global Cafeed epsilonbed k qsaturation L ros q
Cafeed=0.014158;% concentration [mol/dm3]
epsilonbed=0.4; % bed void [-]
portata=102.6522/60;% feed flow rate [dm3/min]
A=3.14*0.238^2/4; %bed area [dm2]
v=portata/A; % velocity [dm/min]
k=0.000156; % costant [1/min]
qsaturation=7.293467/44; % saturation condition [mol CO2/g zeolite]
L=1.1245; % led leight [dm]
ros=2000; % density [g/dm3]
nx = 100; % to define the dimension of the initial condition
y0 = zeros(nx,1); % Initialization of boundary condition
y0(1) = Cafeed;% Boundary condition
xstart = 0.0;% To define the geometry of the axial direction
xend = L;
x = linspace(xstart,xend,nx);
tstart = 0.0; % To define the time range
tend = 30;
tspan = linspace(tstart,tend,100);
[t,y] = ode15s(@(t,y)fun(t,y,x,v),tspan,y0);
function dy=fun(t,y,x,v)
global Cafeed epsilonbed qsaturation L ros k q
nx = numel(y);% dy has the size equal to y
dy = zeros(nx,1);% Initialization of dy
y(1,1)=Cafeed;
for i = 2:nx % To solve the space dependence using a finite difference
dy(i,1) = -v*(y(i)-y(i-1))/(x(i)-x(i-1))-((1-epsilonbed)/epsilonbed)*ros*k*(qsaturation-fun1(t));
end
end
function q = fun1( t )
syms q(t)% To define the function
k=0.000156;% costant [1/min]
qsaturation=7.293467/44;% From experimental data [mol CO2/g zeolite]
ode=diff(q,t)-k*(qsaturation-q)==0;% To define the ODE
cond=q(0)==0;% Initial condition
qsol(t)=dsolve(ode,cond)%To solve the differential equation [mol CO2/g zeolite]
figure
ezplot(qsol(t),[0,60])% To draw the obtained equation
xlabel('Time [min]');
ylabel('qCO2 [mol/g]');
end

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

回答（1 个）

Walter Roberson
Walter Roberson 2018-11-30
function dy=fun(t,y,x,v)
That is going to bring a numeric t into fun.
for i = 2:nx % To solve the space dependence using a finite difference
dy(i,1) = -v*(y(i)-y(i-1))/(x(i)-x(i-1))-((1-epsilonbed)/epsilonbed)*ros*k*(qsaturation-fun1(t));
end
The numeric t is going to get passed into fun1
function q = fun1( t )
The numeric t is going to get received in fun1
syms q(t)% To define the function
This is equivalent to
t = sym('t');
q = symfun( sym('q(t)'), t);
so after that step, t will have been replaced inside the function with a symbolic t, so the input to the function will have been ignored and q will have been assigned a symbolic function.
k=0.000156;% costant [1/min]
qsaturation=7.293467/44;% From experimental data [mol CO2/g zeolite]
ode=diff(q,t)-k*(qsaturation-q)==0;% To define the ODE
cond=q(0)==0;% Initial condition
qsol(t)=dsolve(ode,cond)%To solve the differential equation [mol CO2/g zeolite]
figure
ezplot(qsol(t),[0,60])% To draw the obtained equation
xlabel('Time [min]');
ylabel('qCO2 [mol/g]');
So the same dsolve() is executed each time and plotted, leaving q unchanged.
So the symbolic q(t) function will be returned out of fun1 into the dy(i,1) expression, so the right hand side of the assignment is going to be a symbolic result. But
dy = zeros(nx,1);% Initialization of dy
it is being put into a numeric location. That is going to fail.
  2 个评论
bianchi fiammetta
bianchi fiammetta 2018-11-30
I have understood my script bug, but how can I solve the error? I should convert from symbolic function to a numeric one. If yes, which command should I use?
thank you
Walter Roberson
Walter Roberson 2018-11-30
We do not know what computation fun1 is intended to solve. Should the passed in t value be used as the initial condition for q(0) ? Should the passed in t value be used for the initial condition, q(passed_t) = 0 ? Should the ode be solved as-is and then evaluated at the passed-in t value? If so then you should only do the dsolve() once, getting out a formula, and thereafter you would evaluate the formula at t; you might want to matlabFunction() the result of the dsolve()

请先登录，再进行评论。

类别

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

标签

Community Treasure Hunt

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

Start Hunting!

Translated by