画出来的图怎么变成直线和不见了？？？
在模型中IKK 从0.09-9变化oed_a 中储存了对应变化的数值， oed_a中各个数值不等，我用断点检测了前十个 结果大致如下
1.0e+021*（2.68 4.48 3.80 3.88 0.21 0.22 0.52 1.69 0.17 0.15）
于是命令plot(IKK, oed_a, r-)
可为什么画出来的线是直线 或直接看不见，oed_d, oed_e,oed_me情况一样 是不是编程问题？？如果是， 应该怎么改！ 谢谢您的帮助！
具体编程如下：clear all
clc;
format long;
IKK_normal=0.9;% initional value of N
for KN=1:100
IKK=IKK_normal*KN/10;
Xinit=';
k=';
sys = mdlInit(@ODESenfun, k , Xinit, zeros(11,10));
Tspan = 0:60:6000;
smltn = smltnInit(@ode15s, Tspan);
n = length(Xinit);
m = length(k);
= mdlPlay(sys, smltn);
S2 = S2DParam(S3);
FIM=S2'*S2;
=eig(FIM);
for j=1:11
lambda(j)=lam(j,j);
end
%############################optimal##############
oed_a(KN)=trace(inv(FIM));
oed_d(KN)=det(FIM);
oed_e(KN)=min(lambda);
oed_me(KN)=max(lambda)/min(lambda);
%########PLOTING########
end
figure
plot(IKK,oed_a,'r-');
hold on;
=min(oed_a);
IKK1=b/10*IKK_normal*ones(size(oed_a));
plot(IKK1,oed_a,'b--');
text(0.1,2e-3,['IKK=' num2str(b/10*IKK_normal) '\muM']);
xlabel('Initial concentration of IKK (\muM)');
ylabel('Trace(FIM^-^1)');
legend('A-optimal criterion');
figure
plot(IKK,oed_d,'r-');
hold on;
=max(oed_d);
IKK1=b/10*IKK_normal*ones(size(oed_d));
plot(IKK1,oed_d,'b--');
text(0.1,10e25,['IKK=' num2str(b/10*IKK_normal) '\muM']);
xlabel('Initial concentration of IKK (\muM)');
ylabel('Det(FIM)');
legend('D-optimal criterion');
figure
plot(IKK,oed_e,'r-');
hold on;
=max(oed_e);
IKK1=b/10*IKK_normal*ones(size(oed_e));
plot(IKK1,oed_e,'b--');
text(0.1,5e6,['IKK=' num2str(b/10*IKK_normal) '\muM']);
xlabel('Initial concentration of IKK (\muM)');
ylabel('\lambda_m_i_n(FIM)');
legend('E-optimal criterion');
figure
plot(IKK,oed_me,'r-');
hold on;
=min(oed_me);
IKK1=b/10*IKK_normal*ones(size(oed_me));
plot(IKK1,oed_me,'b--');
text(0.1,8e5,['IKK=' num2str(b/10*IKK_normal) '\muM']);
xlabel('Initial concentration of IKK (\muM)');
ylabel('\lambda_m_a_x(FIM)/\lambda_m_i_n(FIM)');
legend('modified E-optimal criterion');
相关子程序：（以下保证正确）function mdl = mdlInit(ode, theta, x_init, S_init)
% MDLINIT Initialise the parameters used in the model
% Validate input/output arguments
if nargin~=4
error('Function requires 4 input arguments');
end
% Initialise model data structure
mdl.ode = ode;
mdl.theta = theta;
mdl.x_init = x_init;
mdl.S_init = S_init;
mdl.thetaAeq = [];
mdl.thetaBeq = [];
function = mdlPlay(mdl, smltn)
% MDLPLAY Simulate the model over the time span
%
% Returns the simulated states and the sensitivity matrix
%
% = mdlPlay(mdl, smltn)
% mdl model data structure
% smltn simulation data structure
% X 2D matrix i,j (states, time) NB this is the transpose of normal
% return from the solver
% S 3D matrix i,j,k (states, parameters, time)
% Validate the input/output arguments
if nargin~=2,
error('Two input arguments are required');
end
if nargout<2 || nargout>3
error('Two or three output arguments are required');
end
% Initialise the parameters
num_states = length(mdl.x_init);
num_param = length(mdl.theta);
% Reshape the initial parameter/sensitivity values for the ODE simulation)
xs_init = ;
% Run the ODE simulation
= feval(smltn.solver, mdl.ode, smltn.tspan, xs_init, [], mdl.theta);
num_time = length(t);
% Reshape the state matrix (states, time)
X = X_S(:,1:num_states)';
% Reshape the sensitivity derivatives (states, parameters, time)
if nargout==3,
S = reshape(X_S(:,num_states+1:end)', num_states, num_param, num_time);
end
function dx = ODESenfun(t, x, k)
% ODESENFUN Solve both system dynamic ODEs and Sensitivity Equations
% System Dynamic ODEs
% x_dot = F*theta
% S_dot = J*S + F
% Check input and output arguments are correctly supplied
if nargin~=3
error('Three input arguments are required');
end
numStates = 10;
numPara = 11;
% Calculate the state derivatives
% Need to investigate a way to automatically generate these quantities (from a symbolic model)
%dx(1)=-k1*x(1)*x(6)+km1*x(2)+k4*x(4)-km4*x(1)*x(9)+k6*x(5);
%dx(2)=k1*x(1)*x(6) - km1*x(2) - k2*x(2) + km2*x(3)*x(7);
%dx(3)=k2*x(2) - km2*x(3)*x(7) - k3*x(3)*x(8) + km3*x(4) - k5W*x(3) + km5*x(5);
%dx(4)=k3*x(3)*x(8) - km3*x(4) - k4*x(4) + km4*x(1)*x(9);
%dx(5)=k5W*x(3) - km5*x(5) - k6*x(5);
%dx(6)=-k1*x(1)*x(6) + km1*x(2);
%dx(7)=k2*x(2) - km2*x(3)*x(7);
%dx(8)=-k3*x(3)*x(8) + km3*x(4);
%dx(9)=k4*x(4) - km4*x(1)*x(9);
%dx(10)=k6*x(5);
F = [-x(1)*x(6) x(2) 0 0 0 0 x(4) -x(1)*x(9) 0 0 x(5);
x(1)*x(6) -x(2) -x(2) x(3)*x(7) 0 0 0 0 0 0 0;
0 0 x(2) -x(3)*x(7) -x(3)*x(8) x(4) 0 0 -x(3) x(5) 0;
0 0 0 0 x(3)*x(8) -x(4) -x(4) x(1)*x(9) 0 0 0;
0 0 0 0 0 0 0 0 x(3) -x(5) -x(5);
-x(1)*x(6) x(2) 0 0 0 0 0 0 0 0 0;
0 0 x(2) -x(3)*x(7) 0 0 0 0 0 0 0;
0 0 0 0 -x(3)*x(8) x(4) 0 0 0 0 0;
0 0 0 0 0 0 x(4) -x(1)*x(9) 0 0 0;
0 0 0 0 0 0 0 0 0 0 x(5);];
dx = F*k;
% Calculate the sensitivity derivatives
% Need to investigate a way to automatically differentiate these quantities (from a symbolic model)
S = reshape(x(numStates+1:end), numStates, numPara);
J = [-k(1)*x(6)-k(8)*x(9) k(2) 0 k(7) k(11) -k(1)*x(1) 0 0 -k(8)*x(1) 0;
k(1)*x(6) -k(2)-k(3) k(4)*x(7) 0 0 k(1)*x(1) k(4)*x(3) 0 0 0;
0 k(3) -k(4)*x(7)-k(5)*x(8)-k(9) k(6) k(10) 0 -k(4)*x(3) -k(5)*x(3) 0 0;
k(8)*x(9) 0 k(5)*x(8) -k(6)-k(7) 0 0 0 k(5)*x(3)k(8)*x(1) 0;
0 0 k(9) 0 -k(10)-k(11) 0 0 0 0 0;
-k(1)*x(6) k(2) 0 0 0 -k(1)*x(1) 0 0 0 0;
0 k(3) -k(4)*x(7) 0 0 0 -k(4)*x(3) 0 0 0;
0 0 -k(5)*x(8) k(6) 0 0 0 -k(5)*x(3) 0 0;
-k(8)*x(9) 0 0 k(7) 0 0 0 0 -k(8)*x(1) 0;
0 0 0 0 k(11) 0 0 0 0 0;];
dS = J*S+F;
dx((numStates+1):(numStates*(numPara+1))) = reshape(dS,numStates*numPara,1);
function S1 = S2DParam(S)
% S2DPARAMReshape the sensitivity matrix, S, to a 2D matrix
% S - 3D matrix (states, param, time)
% S1 - 2D matrix (states*time, param)
S1 = zeros(size(S,1)*size(S,3), size(S,2));
num_states = size(S,1);
num_time = size(S,3);
for i=1:num_states
S1((i-1)*num_time+1:i*num_time, :) = shiftdim(S(i, :, :),1)';
end
function smltn = smltnInit(solver, tspan)
% Initialize the parameters used in the ODE simulation
% Validate the input/output arguments
if nargin~=2
error('2 input arguments must be supplied');
end
% Set the data structure's fields
smltn.solver = solver;
smltn.tspan = tspan;
如果你有时间 不妨帮我看下， 在此我衷心感谢您的帮助与参与！！ 1.水平有限, 2.时间有限
想想若是你, LZ会仔细看吗?
想法简化问题, 可能较容易...
看看 提问的智慧！！！！(发帖前请认真阅读)
http://forum.vibunion.com/forum/viewthread.php?tid=21991
回复 2 # ChaChing 的帖子
谢谢您的帮助 我会仔细学习的 1.总共6个程序, LZ想会有几个人帮忙看!?
2.LZ都知道plot(IKK, oed_a, r-)有问题了, 为何不检查矩阵变量的大小? ode_a为1*100, IKK为1*1, 当然画出直线!
3.个人水平专业有限, 仅能猜测修改, 在plot前加上IKK=IKK_normal*/10; 是会有图的, 不过不知对错!?
4.LZ没发现在求inv时, 一直有warning吗? 建议检查矩阵是否正确?
Warning: Matrix is close to singular or badly scaled.
Results may be inaccurate. RCOND = 1.461129e-020.
5.很多写法实在不够效率, 如lambda=diag(lam); % for j=1:11, lambda(j)=lam(j,j);end
6.看看精华老帖, 一定有所得
页:
[1]