0

問題があります。行列 P( , ) を計算する関数 popmesh1 があります。次に、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 に答える