Search code examples
matlabloopsodenumerical-integration

Matlab ODE45. How to change a function inside it while calling it?

The functions w_d(Y,T) and q(tg,tm) are being used in a function called by the ODE solver described below:

function dPdh=gasstep1(~,P)
global rho_solid W_B S e height_dryer c_a c_b tm_o

dPdh=zeros(size(P));
Y=P(1);tg=P(2);

% Mass Balance
dPdh(1)=rho_solid*(1-e)*S*height_dryer/W_B*w_d(Y,tg);
% Enthalpy balance
dPdh(2)=-rho_solid*(1-e)*S*height_dryer/W_B*(1/(c_b+c_a*Y))...
         *(q(tg,tm_o)+(c_a*(tg-tm_o)*w_d(Y,tg)));

dPdh=[dPdh(1)
      dPdh(2)];

end

The ODE solver is called in the following manner:

tspan=linspace(0,height_dryer,100);
format long
[h,P]=ode45(@(h,P)gasstep1(h,P),tspan,[Y_o tg_o]);

The values from the generated arrays and P are used in a for loop to compute variables:

X1=zeros(1,Nh);Tm=zeros(1,Nh);
X1(1,1)=X_o;Tm(1,1)=tm_o;

for i=1:Nh-1
    % Mass balance
    X1(1,i+1)=X1(1,i)-dt*((P(i+1,1)-P(i,1))/(h(i+1)-h(i)))*(W_B/S)*(1/(rho_solid*(1-e)));
    % Enthalpy balance
    Tm(1,i+1)=Tm(1,i)+dt*(1/(c_s+c_al*X1(1,i)))*(q(((P(i+1,2)+P(i,2))/2),Tm(1,i))-...
              ((c_a-c_al)*(Tm(1,i))+hv(Tm(1,i)))*w_d((P(i+1,1)-P(i,1)),     ((P(i+1,2)+P(i,2))/2)));
end

The doubt I have is whether the variables in the function w_d and q also change when I call the ODE solver in a loop. I found a method to implement the ODE solver in a loop on the mathworks website (link). I have found help provided when ODE parameters need to be changed (link).

I think my final code should look something like:

for j=2:Nt
    for i=1:Nh-1
        % compute X1 & Tm
    end
    % change variables in function w_d and q
    % [call the ODE]
end

Could anyone please give me an idea on how to start the loop which would help me to complete the code.

Solution

  • The links you pointed to should answer your question: you could make w_d and q parameters of the gasstep function, and use the function parameterization as in the second link you posted.

    From the code you posteed, I understand that you want w_d and q to be dependent on X1 and Tm, right? I'm going to assume that's the case. If they depend on other parameters, just use those instead of X1 and Tm in my explanation:
    make gasstep1 dependent on X1 and Tm and pass them along to w_d and q: your gasstep1 function now has 4 arguments: gasstep1(h, P, X1, Tm), and w_d and q also (w_d(Y, tg, X1, Tm) for instance). Then you only need to create gasstep2 as such:

    gasstep2 = @(h,P) gasstep1(h,P,X1,Tm)

    and you can call ode45 on gasstep2!

    Be careful that in your final code, if you write it like you presented above, X1 and Tm are going to depend on P, which needs to be initialized!