How do I solve Time dependent parameter in ODE

Hello, I am facing problem in solving time dependent parameter in ODE solver. In my dydt fuction, I want G parameter to change as t changes in dydt equation. How can i give input of G parameter in dydt equation. I have tried the below code but it shows error.

[t, ymodel] = ode45(@DiffEqs_crytallization, myODE, tspan, y0, options)

where myODE: (G is of the same matrix size as t and I want to use G in my main function @DiffEqs_crystallization))

function [G] = myODE(t, S)
S=[1.0032; 1.1425; 1.323; 1.53; 1.28; 1.22; 1.09; 1.08];
Kg=exp(3.82);
g=1.62;
G=Kg.*((S-1).^g);


function [dydt] = DiffEqs_crytallization(t,y)
for i= 1:N
       for j=2:N
              extravector1=extravector1+ (y(i)./xmd(j));
       end
         %dydt(i)=(G./delx(i))*(0-0.5*(y(i+1)+y(i)))+k*(xmd(i+1)-xmd(i))*extravector1;
end
dydt=dydt';
return;

[Stephan]

How is G time dependent here? The way it is coded makes it a vector of constant values. No time dependency.

[Jyoti Yadav]

Hello Stephan, thanks for your response.

yes G is a vector of constant value of size same as the time span.

e.g ,, for tspan=[0:10:70]

at t=0, G has some constant value then t=10 G has some another constant value and so on.

but when I use G in dyty equation, it gives error.

function [dydt] = DiffEqs_crytallization(t,y)

for i= 1:N

for j=2:N

extravector1=extravector1+ (y(i)./xmd(j));

end

dydt(i)=(G./delx(i))*(0-0.5*(y(i+1)+y(i)))+k*(xmd(i+1)-xmd(i))*extravector1;

end

dydt=dydt';

return;

error:

Unable to perform assignment because the left and right sides have a different number of elements.

Answer 1

[Alan Stevens]

Perhaps your code needs to be structured more along the following lines;

tspan = .....
y0 = .....
options = .....

[t, ymodel] = ode45(@myODE, tspan, y0, options);

function dydt = myODE(t, y)
    S=[1.0032; 1.1425; 1.323; 1.53; 1.28; 1.22; 1.09; 1.08];
    Kg=exp(3.82);
    g=1.62;
    G=Kg.*((S-1).^g);  % If G is a function of t then express that here
    delx = ... % Needs to be defined
    xmd = .... %  Ditto
    k = ...    % Ditto
    extravector1 = ...% Needs to be initialised
    
    for i= 1:N
       for j=2:N
              extravector1=extravector1+ (y(i)./xmd(j));
       end
       dydt(i)=(G./delx(i))*(0-0.5*(y(i+1)+y(i)))+k*(xmd(i+1)-xmd(i))*extravector1;
    end
        dydt=dydt';

end

[Jyoti Yadav]

Yes Alan, I have included all the values.

global x % particle size, upper edge of each size class, compatible with cumulative passing formulation

global xmax % maximum size of the partcile

global N % number of size classes

global G

global Kg % growth rate coefficient

global g % growth order

global Kj1 % Primary nuleation parameter

global Kj2 % Primary nuleation parameter

global Kb1 % induced nuleation rate coefficient

global b1 % Induced nuleation rate order

global J

global k % breakage rate

global b % breakage distribution function

global delxi

global delx

global xmd

S=[1.0032; 1.1425; 1.323; 1.53; 1.28; 1.22; 1.09; 1.08];

C=[0.079; 0.079; 0.079; 0.079; 0.056; 0.043; 0.035; 0.031];

Cs=[0.0721; 0.063; 0.055; 0.048; 0.042; 0.036; 0.031; 0.028];

Cs3=Cs.*1000000/150.148;

xmax=1000;

x0=1;

N=30;

R=1.26;

% A geometric progression of R

x(1)=1;

for i =2:N

x(i)=R*x(i-1); % geometric progression

end

x(N+1)=1000;

%Growth term

Kg=exp(3.82);

g=1.62;

G=Kg.*((S-1).^g)

%Nuleation

Kj1=exp(24.62);

Kj2=0.0204;

Kb1=exp(24.97);

b1=1.22;

Cc=10442;

Jprim=[0; 2.24e-6; 3894488; 6.37e+8; 7494.39; 1.3e-25; 4.88e-65; 3.77e-91];

%Jprim=Kj1.*S.*exp((-Kj2*(((log(Cc./Cs3)).^3)./((log(S)).^2))));

Jind=Kb1.*((S-1).^b1);

J=Jind+Jprim

% Prediction horizon:

t0=0.0; % Initial time t0

tf=70.0; % Final simulation time

tint=10.0; % Print interval

tspan=[t0:tint:tf]';

options=odeset('RelTol',1e-4,'AbsTol',1e-6,'MaxStep',0.1);

meansize=100.0;

standarddev=50.0;

%% generate the size vector with mid-size of each size class

for i = 1:N

xmd(i) = (x(i)*x(i+1))^0.5; % geometric

end

%

%% generate the density function

density = normpdf(xmd,meansize,standarddev);

for i=1:N

delx(i)=x(i+1)-x(i);

end

%% convert the density to frequency

y0approximate=density.*delx;

sum1=sum(y0approximate);

%Normalize to make the total mass fraction equal 1.0

y0=y0approximate/sum1;

plot(xmd,y0)

%% Call the ODE solver:

[t,ymodel] = ode15s(@DiffEqs_crytallization,tspan,y0,options)

function [dydt] = DiffEqs_crytallization(t,y)

% Declare the global variables

global x % particle size, upper edge of each size class, compatible with cumulative passing formulation

global xmax % maximum size of the partcile

global N % number of size classes

global G

global Kg % growth rate coefficient

global g % growth order

global Kj1 % Primary nuleation parameter

global Kj2 % Primary nuleation parameter

global Kb1 % induced nuleation rate coefficient

global b1 % Induced nuleation rate order

global J

global k % breakage rate

global b % breakage distribution function

k=0.03;

global delx

global xmd

e=10^-10;

extravector1=0;

for i= 1:N

if i==1

for j=2:N

extravector1=extravector1+ (y(i)./xmd(j));

end

dydt(i)=(J./delx(i))+(G/delx(i))*(0-0.5*(y(i+1)+y(i)))+k*(xmd(i+1)-xmd(i))*extravector1;

elseif i<N

extravector2=0.0;

for j= i+1:N

extravector2 = extravector2+y(j)/xmd(j);

end

r(i)=(y(i+1)-y(i)+e)/(y(i)-y(i-1)+e);

phi(i)=max(0,min(2*r(i),min((1/3)+((2*r(i))/3),2)));

dydt(i)=(G/delx(i))*(y(i-1)+0.5*phi(i)*(y(i-1)-y(i-1))-(y(i)+0.5*phi(i)*(y(i)-y(i-1))))+k*(x(i)-x(i-1))*(y(i)/x(i))+k*(x(i+1)-x(i-1))*extravector2-k*y(i);

else

r(i)=(y(i)-y(i)+e)/(y(i)-y(i-1)+e);

phi(i)=max(0,min(2*r(i),min((1/3)+((2*r(i))/3),2)));

dydt(i)=(G/delx(i))*(y(i-1)+0.5*r(i)*(y(i-1)-y(i-2))-(y(i)+0.5*(y(i)-y(i-1))))+k*(x(i)-x(i-1))*(y(i)/x(i))-k*(y(i));

end

end

dydt=dydt';

return;