changing number of variables in ODE

2 views (last 30 days)
H.O.Z.
H.O.Z. on 27 Jun 2015
Commented: Walter Roberson on 27 Jun 2015
Hi All,
I need to write a function that describes the diffusion of materials in a two-layer sphere (a shell layer and a core layer) into the gas. I use the ode functions in matlab to model the concentrations in gas, shell, and core.
My main function is:
r01 = 8; % initial thickness of the core
r02 = 2; % initial thickness of the shell
d0 = (r01+r02)*2; % initial diameter of the sphere
SA0 = zeros(2,1);
V0 = zeros(2,1);
SA0(1) = pi*d0^2; % initial surface area of the sphere
SA0(2) = pi*(2*r01)^2; % initial surface area the core
V0(1) = pi/6*(d0^3-(2*r01)^3); % initial volume of the shell
V0(2) = pi/6*(2*r01)^3; % initial volume of the core
c0 = [0;V0(1)*1e5;V0(2)*1e5];% initial concentrations
options = odeset('RelTol',1e-10);
[t,c] = ode15s(@myfun,[0:120],c0,options,V0,c0);
subplot(1,2,2)
plot(t,c(:,1),'r.',t,c(:,2),'b.',t,c(:,3),'g.');
The function code is:
function y = myfun(t,c,V0,c0)
V=zeros(2,1);
r=zeros(2,1);
SA=zeros(2,1);
V(1) = V0(1)*c(2)/c0(2);% new volume of the shell
V(2) = V0(2)*c(3)/c0(3);% new volume of the core
d = (6*sum(V)/pi)^(1/3);% new diameter of the sphere
r(2) = (6*V(2)/pi)^(1/3)/2;% new thickness of the core
r(1) = d/2 - r(2);% new thickness of the shell
SA(1) = pi*d^2;% new surface area of the sphere
SA(2) = pi*(2*r(1))^2;% new surface area of the core
y(1) = 2e-6*(1e9-c(1))*SA(1); % gas
y(2) = -2e-6*(1e9-c(1))*SA(1); % shell
y(3) = 0; % core
y = y(:);
subplot(1,2,1)
plot(t, r(1),'ro');hold on;
The thickness of the shell get less than zero at some point, which is not reasonable. I want to force the two layers mix into one layer when the shell thickness reaches below 0.2. Then the materials in the new single layer diffuse into gas at the same rate as the initial shell layer. Any idea how to do that? I understand that this modification will change the number of variables from the initial condition (3 variables into 2 variables). Does ode15s allow this? Thank you so much.
  2 Comments
Walter Roberson
Walter Roberson on 27 Jun 2015
Please rewrite using @(t,y) myfun(t, y, v0,c0) instead of @myfun,and remove v0 and c0 from the ode15 calling parameters. There are some odd situations that can crop up with the old form that you use.
Walter Roberson
Walter Roberson on 27 Jun 2015
You should use a plot function in your options instead of plotting in your myfun

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!