How to use nlinfit for a function with a nested infinite sum?

3 次查看（过去 30 天）
Ishan Kulkarni
Ishan Kulkarni 2016-5-15
评论： Clay Swackhamer 2020-8-28
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:
clear
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);
else
S(a) = (1/((2*(n)+1)^2))*exp(((-D)*((2*(n)+1)^2)*(pi^2)*(t))/L^2);
end
end
M(t) = 1-((8/pi^2)*S(20)); %// This is the mass fraction (the y-value)
clear S
% disp(t);
end
figure;
plot(M);
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:
clear
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,x) 1-((8/pi^2)*((1/((2*(n)+1)^2))*exp(((-b(1)).*((2*(n)+1)^2)*(pi^2)*(x))./L^2)));
% 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.
beta0= 1.2*10^-15;
end
[beta,R,J,CovB,MSE]= nlinfit(x,y,modelfun,beta0);
figure
plot(x,y,'s','markersize',5,'color',[0,0,0]);
hold on
xf = linspace(x(1), x(length(x)));
plot(xf,modelfun(beta,xf),'linewidth',1,'color',[1,0,0]);
legend('original data','fit data','location','Best') % the result
end
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.

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

回答（2 个）

Chris J
Chris J 2020-3-2
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?
Clay Swackhamer
Clay Swackhamer 2020-8-28
Hi Ishan and Chris,
I attached an example. If you figured this out already I would be very grateful to see how you did it, since I use this approach in my research as well. Disciples of Peppas unite!

类别

在 Help Center 和 File Exchange 中查找有关 Particle & Nuclear Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by