How to use nlinfit for a function with a nested infinite sum?
6 views (last 30 days)
For my research, I am trying to model (and fit data to) the drug release profile for a drug through a hydrogel, to find the hydrogel diffusivity constant (D). The analytical formula given by the book is:
which, from my understanding, shows the fraction of drug that has released into the environment, and iterates through the infinite sum for each time (t). This yields a graph in the form:
I am supposed to fit experimental data to find D.
Initially I tried to see if I could generate such a graph by just plugging in a value for D using the following code:
D = 1.2*10^-15; %// Assigns value for Diffusivity (D); units=(cm^2/s)
L = 1/10000; %// Assigns value for Thickness of Slab (L); units=cm
for t = 1:1209600 %// defines a time range; units= (s)
S = zeros(20,1); %// defines summation below for set number of terms
for a = 1:200
n = a-1;
if a > 1
S(a) = S(n) + (1/((2*(n)+1)^2))*exp(((-D)*((2*(n)+1)^2)*(pi^2)*(t))/L^2);
S(a) = (1/((2*(n)+1)^2))*exp(((-D)*((2*(n)+1)^2)*(pi^2)*(t))/L^2);
M(t) = 1-((8/pi^2)*S(20)); %// This is the mass fraction (the y-value)
which generates a figure that seems to be correct:
The issue I run into is when I try to use nlinfit to fit experimental data to this equation. I realize that I probably need to define an anonymous function which I can pass on to nlinfit, but I am unsure of how to do that, as it seems I need to define a function within a function. All the documentation I have read deals with functions that do NOT have something that needs to be summed (iterated) for each x (time value).
The diffusivity (D) is unknown, and the coefficient that I'm solving for. I initially tried to define the S function from the earlier code as a function Z(t), which I could then pass onto nlinfit, but it returns an error because D does not have a known value. I reverted my code back to something that at least gave me an output, even though I know it is wrong. Here is what my code looks like:
cd(fileparts(mfilename('fullpath'))) %// Changes the directory to the current location of the script. This makes it easy to just make a folder with script & place excel.
x= xlsread('ReleaseData.xlsx', 'Sheet1', 'A2:A11'); %// reads time data (t) from Excel Sheet
y= xlsread('ReleaseData.xlsx', 'Sheet1', 'B2:B11'); %// reads M/M data (y) from Excel Sheet
%// Here is the data that is being pulled from the sheet:
% x y
% 2 0.081102667
% 4 0.596991
% 8 4.666506
% 24 14.83437333
% 48 23.223108
% 72 29.575486
% 96 32.873329
% 120 33.08912267
% 144 33.24061033
% 168 33.49927467
L = 1/10000;
S = zeros(20,1);
for a = 1:200
n = a-1;
% modelfun=@(b,t) Z(t); // fragment of deleted code where I was trying Z(t) as a previously defined anonymous fxn, but had no value that I could use for D.
xf = linspace(x(1), x(length(x)));
legend('original data','fit data','location','Best') % the result
This gives a beta0 of 1.17e-09 which is many orders of magnitude off from the projected (e-15) which is around where the diffusivity for drugs through hydrogel should be, and causes the analytical function defined in the first part of the post, if given that value for D to instantly go from 0 to 1. The output for the graph from this most recent code gives this:
which does nothing more than show that the fit line is flat (single value) and that the code is wrong.
I apologize if this is long winded but any help or input is greatly appreciated as my research is completely held up by my inability to fit the data to the model equation.
Thank you very much in advance.
Chris J on 2 Mar 2020
Hi Ishan, I'm dealing with the similar equation. Using nlinfit get the same result as your result. Once i use symsum, this kind of problems happened. Have you solved this problem or had any clue about this?