How do I solve Time dependent parameter in ODE

2020 10 21
1 个回答
回答已采纳
更新时间：2020 10 22
1 次查看（30 天）

0 个投票

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;

2 个评论
显示 无 隐藏 无

Stephan
Stephan 2020-10-21
编辑：Stephan 2020-10-21
How is G time dependent here? The way it is coded makes it a vector of constant values. No time dependency.
Jyoti Yadav
Jyoti Yadav 2020-10-22
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.

请先登录，再进行评论。

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

 采纳的回答

Alan Stevens
Alan Stevens 2020-10-21
编辑：Alan Stevens 2020-10-21

0 个投票

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

1 个评论
显示 -1更早的评论 隐藏 -1更早的评论

Jyoti Yadav
Jyoti Yadav 2020-10-22
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;

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by