How to solve differential equations for variable input parameters using for loop and ode45?

19 views (last 30 days)
Matthew Mergy
Matthew Mergy on 12 Aug 2020
Answered: Ayush Gupta on 3 Sep 2020
I have written code (shown below) that uses ode45 to solve 2 differential equations to determine particle number during perikinetic and orthokinetic aggregation for a given flow rate (Reynolds number) through a static mixer. The final particle count from the first equation serves as the initial condition for the second equation.
% Length
L_1 = linspace(0,Trans_L,1000); % cm
L_2 = linspace(L_1(end,1000),L,1000); % cm
% Time
t_1 = (L_1./v)'; % s
t_2 = (L_2./v)'; % s
% Perikinetic Section
N0_1 = [N0; 0];
[t_1,N_1] = ode45(@(t,N) PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G),t_1,N0_1);
% Orthokinetic Section
N0_2 = [0,N_1(end,1)];
[t_2,N_2] = ode45(@(t,N) PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G),t_2,N0_2);
function dNdt = PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G)
dNdt = zeros(2,1);
dNdt(1) = -KA * N(1)^2;
dNdt(2) = -(4/pi) * Ortho_alpha * Vol_frac * G * N(2);
end
Now I am trying to write code that will solve these differential equations for multiple different Reynolds numbers (Re). I have the following code to determine time for each section:
% Length
for j = 1:1:length(Re)
if Trans_L(j) > L
L_1(j,:) = linspace(0,L,1000);
else
L_1(j,:) = linspace(0,Trans_L(j),1000);
end
L2_(j,:) = linspace(Trans_L(j),L,1000);
end
% Time
for j = 1:1:length(Re)
t_1(:,j) = (L_1(j,:)./v(j)); % s
t_2(:,j) = (L_2(j,:)./v(j)); % s
end
But I am unsure how to proceed from here. Should a for loop be used in the function or in the main code? How is indexing properly utilized?
Note: For the inputs to the function, G will vary with different Reynolds number values while the other inputs (KA, Ortho_alpha, Vol_frac) are constants.
Thanks!

Answers (1)

Ayush Gupta
Ayush Gupta on 3 Sep 2020
The code that will solve these differential equations for multiple different Reynolds numbers (Re) can be written as follows and the above script can be broken down into two functions mentionded afterwards.
for j = 1:1:length(Re)
if Trans_L(j) > L
L_1(j,:) = linspace(0,L,1000);
else
L_1(j,:) = linspace(0,Trans_L(j),1000);
end
L2_(j,:) = linspace(Trans_L(j),L,1000);
t_1(:,j) = (L_1(j,:)./v(j)); % s
t_2(:,j) = (L_2(j,:)./v(j)); % s
%Call the ode solver function for each reynold number
% CAlculate G for each specific Reynold's number since it varies with
% Reyolds number, replace it with how it actually varies
G = Re; % for simplication it is being equated to Re
ode_solver_function(L_1(j,:), L2_(j,:), N0, t,N,KA,Ortho_alpha,Vol_frac,G,);
end
%
%function output = a2(L_1, L_2, v,t,N,KA,Ortho_alpha,Vol_frac,G)
function [t_1,N_1,t_2,N_2] = ode_solver_function(L_1, L_2,N0,KA,Ortho_alpha,Vol_frac,G)
% Length
%L_1 = linspace(0,Trans_L,1000); % cm
% L_2 = linspace(L_1(end,1000),L,1000); % cm
% Time
t_1 = (L_1./v)'; % s
t_2 = (L_2./v)'; % s
% Perikinetic Section
N0_1 = [N0; 0];
[t_1,N_1] = ode45(@(t,N) PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G),t_1,N0_1);
% Orthokinetic Section
N0_2 = [0,N_1(end,1)];
[t_2,N_2] = ode45(@(t,N) PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G),t_2,N0_2);
end
function dNdt = PKM_DE_Initial_Fun(t,N,KA,Ortho_alpha,Vol_frac,G)
dNdt = zeros(2,1);
dNdt(1) = -KA * N(1)^2;
dNdt(2) = -(4/pi) * Ortho_alpha * Vol_frac * G * N(2);
end

Products


Release

R2019a

Community Treasure Hunt

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

Start Hunting!