Function not taking updated value of variable when it is passed to it
    9 views (last 30 days)
  
       Show older comments
    
        figure(1)
        t = tiledlayout(2,1);
        ax1 = nexttile;
        ax2 = nexttile;
        k0 = 0.25;
        tdata = [ 0.25; 1; 2; 4; 8; 12 ];                   % Time points of measured drug [] (in hours)
        Y1 = [ 4954.4; 708.8; 85.6; 15.2; 5.6; 1.8 ];       % Primate 1 data
        Y2 = [ 5176.8; 355.1; 61.2; 28.8; 7.2; 0 ];         % Primate 2 data
        yavg = ((Y1+Y2)./2);                                % Averaged primate data
        ydata = yavg./max(yavg);                            % Normalized averaged primate data
        for i = 1:5
            rng('shuffle');
            kin = 10*rand(1,1);
            % test(i,1) = kin; 
            [kfit, ~, ~, exitflag(i)] = lsqcurvefit(@odesolver,k0,tdata,ydata);
            tfit = linspace(0, max(tdata), 100);
            %{
                If exitflag > 0 : The function converged to a solution x.
                If exitflag = 0 : The maximum number of function evaluations or iterations was exceeded.
                If exitflag < 0 : The function did not converge to a solution.
            %}
            Cfit = odesolver(kfit,tfit);
            % check(:,i) = Cfit;
            s(i) = semilogy(ax1,tfit,Cfit,'-o');
                title(ax1,'1st Order Elimination Model');
                xlabel(ax1,'Time (hours)');
                ylabel(ax1,'Concentration (C)');
            if i == 1, hold(ax1,'on'), end                  
            lgdstr{i} = sprintf('k_e = %f\n', kfit);
            %%% Residual Analysis %%%
            ID = interp1(tfit, 1:length(tfit), tdata, 'nearest');
            for k = 1:length(ID)
                res(k) = ydata(k) - Cfit(ID(k));
                res2(k) = res(k)^2; 
            end
            hold(ax2,'on')
            a(i) = plot(ax2, tdata, res2, '-o');
                title(ax2,'Residuals')
                hline = refline(ax2, 0, 0);
                hline.Color = 'k'; 
            lgdstrRes{i} = sprintf('Residual %d', i);     
        end
            legend(s, lgdstr)
            legend(a, lgdstrRes)
            hold(ax2,'off');
            hold(ax1,'on');
            %%% Plots of Primate Data
            semilogy(ax1,tdata,ydata,'--', 'DisplayName','Primate Data',...
                'LineWidth',2)
            hold off
    end 
    end
    function C = odesolver(kin, tdata)                              
            %test(i,1) = kin; 
            %%% Drug PK in NP %%%
            f = @(tdata,u) (-kin*u);                     % 1st Order Equation
            %%% ODE45 Solver %%%
            ut0 = 1;                                     % Initial compartment concentrations
            [tspan, u] = ode45(f, tdata, ut0);           % Solve 1st Order Eq. for u
            C = u;
    end
Hello, I am trying to use least squared fit to fit my first order model to a test data set. My code uses an ODE45 solver to first solve for the concentration and then passes this values to the least squard fit function which then outputs the best rate constant fit. This fitted rate constant is then plugged back into the ODE45 solver and the output concentration is then plotted versus the test data set. I do this 5 times, however, after the first iteration, the following random values for the intial kin value is not being passed to the function so I end up getting the same kfit from the lsqcurvefit function and Cfit from the ODE45 solver for all 5 iterations. I am not sure why the kin value is not being updated with each iteration. 
0 Comments
Answers (1)
  Ayush Gupta
    
 on 15 Sep 2020
        The problem in the above code is the variable kin that is changing with every iteration, it is not being used anywhere for plotting, which is the main reason why the same graph is being generated.  This can be done while calling odesolver function
Cfit = RP(kin,tfit); 
0 Comments
See Also
Categories
				Find more on Ordinary Differential Equations in Help Center and File Exchange
			
	Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
