0

我遇到了一个问题，我有一个函数 popmesh1 可以计算矩阵 P( , )。然后我使用 sens_analysis 存储该矩阵的元素并继续。如何输出 P( , ) 以便跟踪它？我一直认为矩阵的大小是（0,0）。另外，如何将矩阵传递给另一个函数？抱歉发布整个代码，我想具体而清晰，我对 MATLAB 很陌生

    function pdemeshpop(varargin)
global par;
global re;
global P;
global par_n;
global a1;
clc
clear


%INIITIALIZE MESH:- Can change time and age for refining of mesh
time=linspace(0,800,4000);
age=linspace(0,time(end),4000);

dt=time(2)-time(1);
dtao=dt;

P=zeros(length(time),9); % State matrix over time.

P1=zeros(length(time),length(age)); % Mesh for population of P1mod.


prodrev=zeros(length(time),length(age));
p1tot=zeros(length(time));
p2tot=zeros(length(time));
f=zeros(length(time));
A_1=zeros(length(time),1);


%Parameters
G=log(2)/30; %This growth rate had been set to nthroot(2,20), but I think it should be log(2)/20 for a doubling time of 20 mins. Units 1/min
R=.75; %Reduction in growth rate due to viral production, range from 0.5-0.75
global A_s; %Number of virus produced each minute from one cell? Units 1/min
A_s = 35; %Source for this?
global re;
re = varargin(1); %Reduction in efficiency of virus production in P1mod
c=1.5e9; %Concentration of cells in saturated culture. Units 1/cm^3 Source: http://bionumbers.hms.harvard.edu/bionumber.aspx?&id=100984&ver=2
K=3e-11; %Adsorption rate. Units cm^3/min. Source: Tzagoloff, H., and D. Pratt. 1964. The initial steps in infection with coliphage M13. Virology 24:372?380.
i = 1; %Is flipping of switch induced or not induced? if i==1 then switch is induced.

if i==1
    S_i=0.5; %S_i is probability that a switch will flip during a timestep in a p1ori cell. Units pure number (a probablility). Ranges from 0 to 1.
elseif i==0
    S_i= 0.005;
end

%IC and BC implementation for the 9 dependent variables <<<10?
P0=zeros(9,1);
P0(1)=100; %Initial concentration of senders. Units: cells/ml.
P0(2)=10000; %Initial concentration of primary receivers. Units: cells/ml.
P0(3)=10000;  %Initial concentration of secondary receivers. Units: cells/ml.

%The loop below covers the initial conditions and BC of t=0,all ages
for i=1:9
    P(1,i)=P0(i); 
end


%Iterative solution

for m=1:length(time)-1 % m is timestep

        %Simplifications
        p1tot(m)=P(m,2)+P(m,4)+P(m,6);
        p2tot(m)=P(m,3)+P(m,5)+P(m,7);
        f(m)=1-(P(m,1)+p1tot(m)+p2tot(m))/c; 

        %Senders
        P(m+1,1)=dt*(P(m,1)*G*R*f(m))+P(m,1);

        %Primary Receivers
        P(m+1,2)=dt*((P(m,2)*G*f(m))-K*(P(m,8)+P(m,9))*P(m,2))+P(m,2); 

        %Secondary Receivers
        P(m+1,3)=dt*((P(m,3)*G*f(m))-K*(P(m,8)+P(m,9))*P(m,3))+P(m,3);

        %Primary Original
        P(m+1,4)=dt*((P(m,4)*G*f(m))+K*P(m,8)*P(m,2)-S_i*P(m,4))+P(m,4);

        %Secondary Original
        P(m+1,5)=dt*((P(m,5)*G*f(m))+K*P(m,8)*P(m,3))+P(m,5);

    for n=1:m
        t=(m-1)*dt; %Why not t=m*dt?
        tao=(n-1)*(dtao);%Determines current age basket
        prodrev(m,n)=rate(t-tao); %Calculates corresponding rate of production of phage (reversed)   

        %Primary Modified
        if n==1
        P1(m+1,n)=dt*(K*P(m,2)*P(m,9)+S_i*P(m,4)); %Left hand side boundary (New cells at age zero)  
        else
        P1(m+1,n)=dt*(-((P1(m,n)-P1(m,n-1))/dtao)+P1(m,n)*G*R*f(m))+P1(m,n);    
        end

    end

        P(m+1,6)=sum(P1(m+1,:)); %phi1mod

        %Secondary Modified
        P(m+1,7)=dt*((P(m,7)*G*f(m))+K*P(m,9)*P(m,3))+P(m,7);

        %Original
        P(m+1,8)=dt*(A_s*P(m,1))+P(m,8);



        if m<2
            A_1(m)=0;
        else
            convolution(m,:)=prodrev(m,:).*P1(m,:); 

           %A_1(m)=dtao*trapz(conv(prod1(m,:), P1(m,:))); 
           A_1(m)=dtao*trapz([convolution(m,:) 0]);
           %A_1 obtained by convolving the discrete vectors of P1 and prod1
           %then finding the area under the curve

        end    

        %Modified
        P(m+1,9)=dt*A_1(m)+P(m,9);

end

P1;



end

function prod=rate(tao)
%Function generates production rate values of the infected cells based on their age 
global A_s;
global re;
a=re*A_s; %Max production rate

ageofcell=tao;
if ageofcell<=10
        prod=0;
    elseif ageofcell<=50
        prod=(a/40)*(ageofcell-10);
    else
        prod=a;

end
end

以及调用上述函数的其他代码，并且我想将 P(time(length),7) 传递给：

function sens_analysis

global par;
global re;
global par_n;
global P;
time=linspace(0,800,4000);

pdemeshpop_final_re_sens;
re_0 = re;
par = re;
s_nom_ss = a1;
delta = 0.05;
par_n = par*(1+delta);
pdemeshpop_final_re_sens_par(par_n);  % similar to pdemeshpop_final_re_sens
s_pert_ss = P(length(time),7);
abs_sens = (s_pert_ss - s_nom_ss)/(delta*re_0);
rel_sens = abs_sens*(re_value/s_nom_ss);
end

再次，抱歉发布整个代码，觉得这是一个必要的邪恶。全局变量也可能是不必要的。可能是显而易见的。我可能需要先以某种方式存储 P 。有人可以仔细解释一下吗？谢谢！

4

2 回答 2

0

我在您的代码中注意到的第一件事是您使用了全局变量。如果可以避免，一般不建议这样做。考虑将它们作为函数的输入，单独或在结构中。

当然clear是删除你的变量，但一般来说，如果你想看看发生了什么，试着在你的代码中放置一些断点。这使您可以检查所有现有变量。使用 f10，您可以单步执行代码并查看一切是如何进行的。

此外，我始终建议您使用dbstop if error，这样您可以有效地处理您将遇到的错误。

于 2013-09-19T14:06:52.853 回答
0

在函数的开头，pdemeshpop您有一个clear语句，它正在从您的工作区中删除由全局语句声明的变量。注释掉该clear语句，您将绕过该问题。

于 2013-09-19T08:28:09.237 回答