ODE45 solver returning constant value as the input variable.

2 次查看（过去 30 天）
Pc
Pc 2024-2-16
移动：Torsten 2024-2-16
I'm encountering an issue with my MATLAB code where the first state variable in my system does not seem to evolve over time when solved using an ODE solver. However, the subsequent variables, which are dependent on the first variable, do appear to evolve over time. This inconsistency is puzzling to me. I've tried adjusting the solver tolerances, but the issue persists. I suspect there might be an error in how I've defined the equations or the initial conditions.
Could someone please review my code and help me understand why the first state variable remains constant while the others change over time? Any insights or suggestions would be greatly appreciated. Additionally, please provide explanations in simple language as I'm still learning. Thank you in advance for your assistance!
%Define the variables to be used
%Lc-critical particle length[μm]
%L-characteristic mean particle length[μm]
%Delta_Li= difference between new and previous crystal length [μm]
%Ni-number of particles
%N-sum of N
%n=population density
%h-diffusion layer thickness[μm]
%mo-initial mass of crystals[g]
%mi_j-mass of undissolved crystals[g]
%mw-molecular weight of ibuprofen [g/mol]
%w-particle weight at time t [g]
%KA-surface shape factor
%KV-volume shape factor
%D-diffusion coefficient[cm2/s]
%p- actual crystal density of particles [kg/m3]; will calculate it in g/ml
%Cs-solubility of ibuprofen[g/mL]
%Cb_exp- experimental bulk concentration[g/mL]
%Yexpt-fraction of drug dissolved at time t
%S- surface area[m2/g]
%V= volume of solution[mL]
%v=volume of ibuprofen [mL/mol]
%constants and parameters
mo=1;%[g]
v=200.3;%[mL/mol]
mw=206.28;%[g/mol]
w=1;%[g]
p=mo/((mo/mw)*v);%[g/mL]%calculate density from mass and vol of ibuprofen
p = p * 1e12; % 1e12 is 10^12.% Convert density to g/μm^3
Lc=50;%[μm] %estimated from previous literature according to the article 2
D=8.88*10^-6;%[cm2/s] Averaged from different previous models.
V=100;%[mL] %taken from article 1
KA=pi; %for spherical crystals
KV=pi/6;%for spherical crystals
% Experimental data, Run 6
Cs=0.00425; %[g/g]
% Yexpt mass of drug dissolved at time t/total mass of drug dissolved when dosage form is exhausted
Yexpt= [0.0000 0.1321 0.2447 0.3343 0.4122 0.5298 0.6246 0.6869 0.7348 0.7729 0.7931 0.8114 0.8272]';
tdata=[0 30 60 90 120 180 240 300 360 420 480 540 600]';
%convert Yexpt to Cb_exp
Cb_exp= (Yexpt.*mo)./V;
% initial particle conditions for size 300-500 μm
lowest_L=300;
Upper_L=500;
Lo=(lowest_L+Upper_L)/2; %mean particle size
Delta_L=Upper_L-lowest_L;
n=(mo*w)/(p*KA*(Lo)^3*Delta_L); % population density
Ni=n*Delta_L; % particle number %assumed to remain constant.
mi_j=1; % initial mass of undissolved particles
Cb_cal=0; %at t=0
S=KA * Lo^2 * Ni / V; % intial surface area
dLo=[Lo mi_j Cb_cal S]; % [L mi_j Cb_Cal S]; initial conditions for the ODE solver
h=Lc/2; %constant for a size independent scenario
tspan=tdata;% time over which integration will occur
Function Definitions
% integrate the differential equations
[t, L] = ode45(@(t, L) odefunc(t, L,tdata,Cb_exp, p,KA,h,KV, D, Cs,mo,Ni, V), tspan, dLo);
% Extract Li, mi_j, Cb_cal, and S corresponding to the integrated time points
Li = L(:, 1);
mi_j = L(:, 2);
Cb_cal = L(:, 3);
S = L(:, 4);
% Define the function to calculate SSE
calculate_SSE = @(h) sum_SSE(h, tdata, Cb_exp, p, KA, KV, D, Cs, mo, Ni, V, dLo);
% Perform the calculation of SSE
h_guess=h;
SSE = calculate_SSE(h_guess);
% Perform optimization
options = optimset('Display', 'iter');
h_opt = fminsearch(calculate_SSE, h_guess, options);
disp(['SSE: ', num2str(SSE)]);
disp(['Optimal h: ', num2str(h_opt)]);
figure;
plot(tdata, Cb_exp' ,'r--', 'LineWidth', 2); % Plot Cb_exp in red dashed lines
xlabel('Time (min)');
% Plot Cb_cal and Cb_exp against tdata
figure;
plot(tdata, Cb_cal', 'b-', 'LineWidth', 2); % Plot Cb_cal in blue
hold on;
plot(tdata, Cb_exp' ,'r--', 'LineWidth', 2); % Plot Cb_exp in red dashed lines
xlabel('Time (min)');
ylabel('Concentration (g/mL)');
title('Comparison of Cb\_cal and Cb\_exp');
legend('Cb\_cal', 'Cb\_exp');
grid on;
function dLdt = odefunc(t, L,tdata,Cb_exp, p,KA,h,KV, D, Cs,mo,Ni, V)
% Evaluate experimental concentration at the current time t
[~, index] = min(abs(tdata - t));
Cb_at_t = Cb_exp(index);
%create a place holder for dLdt
dLdt=zeros(4, 1);
%replace dLdt positions with the equations to solve
dLdt(1) = -(KA / (3 * p * KV)) * (D / h) .* (Cs - Cb_at_t); %Equation 1: dL/dt
dLdt(2) = p * KV *(L(1))^3; % Equation 2: mi_j
dLdt(3) = (mo - L(2)) / V; % Equation 3: Cb_cal
dLdt(4) = KA * (L(1))^2 * Ni / V; % Equation 4: S
end
% Define the function to calculate SSE
function sse = sum_SSE(h, tdata, Cb_exp, p, KA, KV, D, Cs, mo, Ni, V,dLo)
[~, L] = ode45(@(t, L) odefunc(t, L, tdata, Cb_exp, p, KA, h, KV, D, Cs, mo, Ni, V), tdata, dLo);
Cb_cal = L(:, 3);
sse = sum((Cb_cal - Cb_exp).^2);
end
  1 个评论
Star Strider
Star Strider 2024-2-16
Ran in:
%Define the variables to be used
%Lc-critical particle length[μm]
%L-characteristic mean particle length[μm]
%Delta_Li= difference between new and previous crystal length [μm]
%Ni-number of particles
%N-sum of N
%n=population density
%h-diffusion layer thickness[μm]
%mo-initial mass of crystals[g]
%mi_j-mass of undissolved crystals[g]
%mw-molecular weight of ibuprofen [g/mol]
%w-particle weight at time t [g]
%KA-surface shape factor
%KV-volume shape factor
%D-diffusion coefficient[cm2/s]
%p- actual crystal density of particles [kg/m3]; will calculate it in g/ml
%Cs-solubility of ibuprofen[g/mL]
%Cb_exp- experimental bulk concentration[g/mL]
%Yexpt-fraction of drug dissolved at time t
%S- surface area[m2/g]
%V= volume of solution[mL]
%v=volume of ibuprofen [mL/mol]
%constants and parameters
mo=1;%[g]
v=200.3;%[mL/mol]
mw=206.28;%[g/mol]
w=1;%[g]
p=mo/((mo/mw)*v);%[g/mL]%calculate density from mass and vol of ibuprofen
p = p * 1e12; % 1e12 is 10^12.% Convert density to g/μm^3
Lc=50;%[μm] %estimated from previous literature according to the article 2
D=8.88*10^-6;%[cm2/s] Averaged from different previous models.
V=100;%[mL] %taken from article 1
KA=pi; %for spherical crystals
KV=pi/6;%for spherical crystals
% Experimental data, Run 6
Cs=0.00425; %[g/g]
% Yexpt mass of drug dissolved at time t/total mass of drug dissolved when dosage form is exhausted
Yexpt= [0.0000 0.1321 0.2447 0.3343 0.4122 0.5298 0.6246 0.6869 0.7348 0.7729 0.7931 0.8114 0.8272]';
tdata=[0 30 60 90 120 180 240 300 360 420 480 540 600]';
%convert Yexpt to Cb_exp
Cb_exp= (Yexpt.*mo)./V;
% initial particle conditions for size 300-500 μm
lowest_L=300;
Upper_L=500;
Lo=(lowest_L+Upper_L)/2; %mean particle size
Delta_L=Upper_L-lowest_L;
n=(mo*w)/(p*KA*(Lo)^3*Delta_L); % population density
Ni=n*Delta_L; % particle number %assumed to remain constant.
mi_j=1; % initial mass of undissolved particles
Cb_cal=0; %at t=0
S=KA * Lo^2 * Ni / V; % intial surface area
dLo=[Lo mi_j Cb_cal S]; % [L mi_j Cb_Cal S]; initial conditions for the ODE solver
h=Lc/2; %constant for a size independent scenario
tspan=tdata;% time over which integration will occur
% integrate the differential equations
[t, L] = ode45(@(t, L) odefunc(t, L,tdata,Cb_exp, p,KA,h,KV, D, Cs,mo,Ni, V), tspan, dLo);
% Extract Li, mi_j, Cb_cal, and S corresponding to the integrated time points
Li = L(:, 1);
mi_j = L(:, 2);
Cb_cal = L(:, 3);
S = L(:, 4);
% Define the function to calculate SSE
calculate_SSE = @(h) sum_SSE(h, tdata, Cb_exp, p, KA, KV, D, Cs, mo, Ni, V, dLo);
% Perform the calculation of SSE
h_guess=h;
SSE = calculate_SSE(h_guess);
% Perform optimization
options = optimset('Display', 'iter');
h_opt = fminsearch(calculate_SSE, h_guess, options);
Iteration Func-count f(x) Procedure 0 1 9.77753e+45 1 2 9.77753e+45 initial simplex 2 5 9.77753e+45 shrink 3 8 9.77753e+45 shrink 4 11 9.77753e+45 shrink 5 14 9.77753e+45 shrink 6 17 9.77753e+45 shrink 7 20 9.77753e+45 shrink 8 23 9.77753e+45 shrink 9 26 9.77753e+45 shrink 10 29 9.77753e+45 shrink 11 32 9.77753e+45 shrink 12 35 9.77753e+45 shrink 13 38 9.77753e+45 shrink 14 41 9.77753e+45 shrink 15 44 9.77753e+45 shrink Optimization terminated: the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-04 and F(X) satisfies the convergence criteria using OPTIONS.TolFun of 1.000000e-04
disp(['SSE: ', num2str(SSE)]);
SSE: 9.777527001913533e+45
disp(['Optimal h: ', num2str(h_opt)]);
Optimal h: 25
figure;
plot(tdata, Cb_exp' ,'r--', 'LineWidth', 2); % Plot Cb_exp in red dashed lines
xlabel('Time (min)');
% Plot Cb_cal and Cb_exp against tdata
figure;
plot(tdata, Cb_cal', 'b-', 'LineWidth', 2); % Plot Cb_cal in blue
hold on;
plot(tdata, Cb_exp' ,'r--', 'LineWidth', 2); % Plot Cb_exp in red dashed lines
xlabel('Time (min)');
ylabel('Concentration (g/mL)');
title('Comparison of Cb\_cal and Cb\_exp');
legend('Cb\_cal', 'Cb\_exp');
grid on;
function dLdt = odefunc(t, L,tdata,Cb_exp, p,KA,h,KV, D, Cs,mo,Ni, V)
% Evaluate experimental concentration at the current time t
[~, index] = min(abs(tdata - t));
Cb_at_t = Cb_exp(index);
%create a place holder for dLdt
dLdt=zeros(4, 1);
%replace dLdt positions with the equations to solve
dLdt(1) = -(KA / (3 * p * KV)) * (D / h) .* (Cs - Cb_at_t); %Equation 1: dL/dt
dLdt(2) = p * KV *(L(1))^3; % Equation 2: mi_j
dLdt(3) = (mo - L(2)) / V; % Equation 3: Cb_cal
dLdt(4) = KA * (L(1))^2 * Ni / V; % Equation 4: S
end
% Define the function to calculate SSE
function sse = sum_SSE(h, tdata, Cb_exp, p, KA, KV, D, Cs, mo, Ni, V,dLo)
[~, L] = ode45(@(t, L) odefunc(t, L, tdata, Cb_exp, p, KA, h, KV, D, Cs, mo, Ni, V), tdata, dLo);
Cb_cal = L(:, 3);
sse = sum((Cb_cal - Cb_exp).^2);
end
.

请先登录，再进行评论。

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

采纳的回答

Torsten
Torsten 2024-2-16
移动：Torsten 2024-2-16
As you can see when you run the below code, p is in the order of 1e12 which makes dL/dt effectivly 0.
Check your units.
%Define the variables to be used
%Lc-critical particle length[μm]
%L-characteristic mean particle length[μm]
%Delta_Li= difference between new and previous crystal length [μm]
%Ni-number of particles
%N-sum of N
%n=population density
%h-diffusion layer thickness[μm]
%mo-initial mass of crystals[g]
%mi_j-mass of undissolved crystals[g]
%mw-molecular weight of ibuprofen [g/mol]
%w-particle weight at time t [g]
%KA-surface shape factor
%KV-volume shape factor
%D-diffusion coefficient[cm2/s]
%p- actual crystal density of particles [kg/m3]; will calculate it in g/ml
%Cs-solubility of ibuprofen[g/mL]
%Cb_exp- experimental bulk concentration[g/mL]
%Yexpt-fraction of drug dissolved at time t
%S- surface area[m2/g]
%V= volume of solution[mL]
%v=volume of ibuprofen [mL/mol]
%constants and parameters
mo=1;%[g]
v=200.3;%[mL/mol]
mw=206.28;%[g/mol]
w=1;%[g]
p=mo/((mo/mw)*v);%[g/mL]%calculate density from mass and vol of ibuprofen
p = p * 1e12; % 1e12 is 10^12.% Convert density to g/μm^3
Lc=50;%[μm] %estimated from previous literature according to the article 2
D=8.88*10^-6;%[cm2/s] Averaged from different previous models.
V=100;%[mL] %taken from article 1
KA=pi; %for spherical crystals
KV=pi/6;%for spherical crystals
% Experimental data, Run 6
Cs=0.00425; %[g/g]
% Yexpt mass of drug dissolved at time t/total mass of drug dissolved when dosage form is exhausted
Yexpt= [0.0000 0.1321 0.2447 0.3343 0.4122 0.5298 0.6246 0.6869 0.7348 0.7729 0.7931 0.8114 0.8272]';
tdata=[0 30 60 90 120 180 240 300 360 420 480 540 600]';
%convert Yexpt to Cb_exp
Cb_exp= (Yexpt.*mo)./V;
% initial particle conditions for size 300-500 μm
lowest_L=300;
Upper_L=500;
Lo=(lowest_L+Upper_L)/2; %mean particle size
Delta_L=Upper_L-lowest_L;
n=(mo*w)/(p*KA*(Lo)^3*Delta_L); % population density
Ni=n*Delta_L; % particle number %assumed to remain constant.
mi_j=1; % initial mass of undissolved particles
Cb_cal=0; %at t=0
S=KA * Lo^2 * Ni / V; % intial surface area
dLo=[Lo mi_j Cb_cal S]; % [L mi_j Cb_Cal S]; initial conditions for the ODE solver
h=Lc/2; %constant for a size independent scenario
tspan=tdata;% time over which integration will occur
% integrate the differential equations
[t, L] = ode45(@(t, L) odefunc(t, L,tdata,Cb_exp, p,KA,h,KV, D, Cs,mo,Ni, V), tspan, dLo);
plot(t,L(:,1))
function dLdt = odefunc(t, L,tdata,Cb_exp, p,KA,h,KV, D, Cs,mo,Ni, V)
% Evaluate experimental concentration at the current time t
%[~, index] = min(abs(tdata - t));
%Cb_at_t = Cb_exp(index);
Cb_at_t = interp1(tdata,Cb_exp,t);
%create a place holder for dLdt
dLdt=zeros(4, 1);
%replace dLdt positions with the equations to solve
dLdt(1) = -(KA / (3 * p * KV)) * (D / h) .* (Cs - Cb_at_t); %Equation 1: dL/dt
dLdt(2) = p * KV *(L(1))^3; % Equation 2: mi_j
dLdt(3) = (mo - L(2)) / V; % Equation 3: Cb_cal
dLdt(4) = KA * (L(1))^2 * Ni / V; % Equation 4: S
end

更多回答（0 个）

类别

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

标签

产品


版本

R2020b

Community Treasure Hunt

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

Start Hunting!

Translated by