Run a system of odes for different parameters and plot
显示 更早的评论
I have a system of 3 odes:
function dydt=odefcnNY_v20(~,y,D,Cs,rho,r0,N,V,k01,F,CL10,V2)
dydt=zeros(3,1);
dydt(1)=(-D*Cs*(1-y(2)))/(rho*r0^2*y(1)); % dr*/dt
dydt(2)=((4*pi*D*N*r0*(1-y(2))*y(1))/V)-(F*k01*y(2)); %dC*/dt
dydt(3)=((k01*y(2)*Cs*V)-(CL10*y(3)))/V2; %dC2/dt
end
I then run this code to simulate the response for a given set of parameters D, Cs, rho, r0, N, V, k01, F, CL10 and V2, to make plots of y(1), y(2), y(3) and various transformations of them and then I compare y(3) to experimental data. Is there a way to run this system of odes in a loop for a different set of k01 values varying from 0.1 to 17000 and to automatically plot the resulting y(3) from the model with the y(3) from the experimental dataset all in one graph, labelled by k01 and clearly labelled by whether it is a model output or an experimental curve?
% Parameters
MW=668.89; % molecular weight
D=9.916e-5*(MW^-0.4569)*60/600000 %m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60], divide by 600,000 to convert to m2/s
rho=1170; %kg/m3
dv90_1M=3.2e-6; %m
r0=dv90_1M/2 %m dv90
Cs=1.18e-3 % kg/m3
V=6.9e-8*60;%m3
W=6.9*1e-6; %kg
N=W/((4/3)*pi*r0^3*rho);
k01=1.7/3600; %s
V2=1.561/1000; % m3
CL10=0.733/3.6e6; % m3/s
F=1;
% initial conditions
tspan=[0 4200*3600]; %
y0=[1 0 0];
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,y]=ode15s(@(t,y) odefcnNY_v20(t,y,D,Cs,rho,r0,N,V,k01,F,CL10,V2), tspan, y0, options);
%Plot y(1) and a transformation of it
plot(t/(3600*24),y(:,1),'-o')%plot time in days, and y(1)
xlabel('time, days')
ylabel('r*, (rp/r0)')
legend('Compound 1')
title ('r*');
plot(t/(3600*24),y(:,1)*r0*1e6); %plot y(1) in microns
xlabel('time, days');
ylabel('r, microns');
legend('Compound 1');
title('r');
%Plot y(2) and a transformation of it
plot(t/(3600*24),y(:,2),'-') %plot time in days, and y(2)
xlabel('time, days')
ylabel('C* (C/Cs)')
legend('Compound 2')
title('C*');
ConcDiss = y(:,2)*Cs; %Concentration dissolved, kg/m3
plot(t/(3600*24), ConcDiss,'DisplayName', 'Theoretical') % time in days, and bulk concentration on y
xlabel('time, days')
ylabel('C1, kg/m3')
title ('C1');
%Plot y(3) and compare to experimental data
plot(t/(3600*24),y(:,3),'-') %plot time in days, and C2
xlabel('time, days')
ylabel('C2')
legend('Compound 1')
title('C2');
hold on
MeasuredData1M=xlsread('PP_Data_1M.xlsx');
t_Measured1M=MeasuredData1M(:,1);
C2_Measured1M=MeasuredData1M(:,2);
plot(t_Measured1M/24,C2_Measured1M,'or','DisplayName','Measured PP1M');
legend('Model', 'Measured PP1M');
2 个评论
Ameer Hamza
2020-3-20
Can you attach the file 'PP_Data_1M.xlsx'? It will help in understanding how you are trying to visualize this data.
Nora Rafael
2020-3-20
采纳的回答
Ameer Hamza
2020-3-20
You can try something like the following code, but there is some issue with your model, which causes ode45 to fail at some point in time.
% Parameters
MW=668.89; % molecular weight
D=9.916e-5*(MW^-0.4569)*60/600000; %m2/s - [D(cm2/min)=9.916e-5*(MW^-0.4569)*60], divide by 600,000 to convert to m2/s
rho=1170; %kg/m3
dv90_1M=3.2e-6; %m
r0=dv90_1M/2; %m dv90
Cs=1.18e-3; % kg/m3
V=6.9e-8*60;%m3
W=6.9*1e-6; %kg
N=W/((4/3)*pi*r0^3*rho);
V2=1.561/1000; % m3
CL10=0.733/3.6e6; % m3/s
F=1;
K01 = linspace(0.1, 17000, 5);
MeasuredData1M=xlsread('PP_Data_1M.xlsx');
t_Measured1M=MeasuredData1M(:,1);
C2_Measured1M=MeasuredData1M(:,2);
% initial conditions
tspan=[0 4200*3600]; %
y0=[1 0 0];
fig = figure();
ax = axes();
hold(ax);
plot(t_Measured1M/24,C2_Measured1M,'or','DisplayName','Measured PP1M');
for i=1:numel(K01)
k01 = K01(i);
options = odeset('RelTol',1e-10,'AbsTol',1e-10);
[t,y]=ode15s(@(t,y) odefcnNY_v20(t,y,D,Cs,rho,r0,N,V,k01,F,CL10,V2), tspan, y0, options);
plot(t/(3600*24),y(:,3),'-', 'DisplayName', ['model k01=' num2str(k01)]); %plot time in days, and C2
end
%Plot y(3) and compare to experimental data
xlabel('time, days')
ylabel('C2')
title('C2');
legend(ax);
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
