画出来的图怎么变成直线和不见了？？？

lphtoby

在模型中
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情况一样 是不是编程问题？？如果是， 应该怎么改！ 谢谢您的帮助！

具体编程如下：

1. clear all
   
2. clc;
   
3. format long;
   
4. 
   
5. IKK_normal=0.9;% initional value of N
   
6. for KN=1:100
   
7.    
   
8.     IKK=IKK_normal*KN/10;
   
9.     Xinit=[1.5e-5 0 0 0 0 0.8 0 IKK 0 0]';
   
10.     k=[100000 1000 100 10000 50000 200 1000 20000 5000 100 500]';
    
11.     sys = mdlInit(@ODESenfun, k , Xinit, zeros(11,10));
    
12.     Tspan = 0:60:6000;
    
13.     smltn = smltnInit(@ode15s, Tspan);
    
14.     n = length(Xinit);
    
15.     m = length(k);
    
16.     [t, X, S3] = mdlPlay(sys, smltn);
    
17. 
    
18.     S2 = S2DParam(S3);
    
19.   
    
20.   FIM=S2'*S2;
    
21. [vec, lam]=eig(FIM);
    
22. for j=1:11
    
23.     lambda(j)=lam(j,j);
    
24. end
    
25. 
    
26. 
    
27. %############################optimal##############
    
28. oed_a(KN)=trace(inv(FIM));
    
29. oed_d(KN)=det(FIM);
    
30. oed_e(KN)=min(lambda);
    
31. oed_me(KN)=max(lambda)/min(lambda);
    
32. %########PLOTING########
    
33. 
    
34.    
    
35. 
    
36. 
    
37. end
    
38. 
    
39. 
    
40. figure
    
41.     plot(IKK,oed_a,'r-');
    
42.     hold on;
    
43.     [a,b]=min(oed_a);
    
44.     IKK1=b/10*IKK_normal*ones(size(oed_a));
    
45.     plot(IKK1,oed_a,'b--');
    
46.     text(0.1,2e-3,['IKK=' num2str(b/10*IKK_normal) '\muM']);
    
47.     xlabel('Initial concentration of IKK (\muM)');
    
48.     ylabel('Trace(FIM^-^1)');
    
49.     legend('A-optimal criterion');
    
50.    
    
51. figure
    
52.     plot(IKK,oed_d,'r-');
    
53.     hold on;
    
54.     [a,b]=max(oed_d);
    
55.     IKK1=b/10*IKK_normal*ones(size(oed_d));
    
56.     plot(IKK1,oed_d,'b--');
    
57.     text(0.1,10e25,['IKK=' num2str(b/10*IKK_normal) '\muM']);
    
58.     xlabel('Initial concentration of IKK (\muM)');
    
59.     ylabel('Det(FIM)');
    
60.     legend('D-optimal criterion');
    
61.    
    
62. figure
    
63.     plot(IKK,oed_e,'r-');
    
64.     hold on;
    
65.     [a,b]=max(oed_e);
    
66.     IKK1=b/10*IKK_normal*ones(size(oed_e));
    
67.     plot(IKK1,oed_e,'b--');
    
68.     text(0.1,5e6,['IKK=' num2str(b/10*IKK_normal) '\muM']);
    
69.     xlabel('Initial concentration of IKK (\muM)');
    
70.     ylabel('\lambda_m_i_n(FIM)');
    
71.     legend('E-optimal criterion');
    
72.    
    
73. figure
    
74.     plot(IKK,oed_me,'r-');
    
75.     hold on;
    
76.     [a,b]=min(oed_me);
    
77.     IKK1=b/10*IKK_normal*ones(size(oed_me));
    
78.     plot(IKK1,oed_me,'b--');
    
79.     text(0.1,8e5,['IKK=' num2str(b/10*IKK_normal) '\muM']);
    
80.     xlabel('Initial concentration of IKK (\muM)');
    
81.     ylabel('\lambda_m_a_x(FIM)/\lambda_m_i_n(FIM)');
    
82.     legend('modified E-optimal criterion');
    
83. 
    

复制代码

相关子程序：（以下保证正确）

1. function mdl = mdlInit(ode, theta, x_init, S_init)
   
2. % MDLINIT   Initialise the parameters used in the model
   
3. 
   
4. % Validate input/output arguments
   
5. if nargin~=4
   
6.     error('Function requires 4 input arguments');
   
7. end
   
8. 
   
9. % Initialise model data structure
   
10. mdl.ode = ode;
    
11. mdl.theta = theta;
    
12. mdl.x_init = x_init;
    
13. mdl.S_init = S_init;
    
14. mdl.thetaAeq = [];
    
15. mdl.thetaBeq = [];
    

复制代码

1. function [t, X, S] = mdlPlay(mdl, smltn)
   
2. % MDLPLAY        Simulate the model over the time span
   
3. %
   
4. % Returns the simulated states and the sensitivity matrix
   
5. %
   
6. % [t, X, S] = mdlPlay(mdl, smltn)
   
7. %   mdl model data structure
   
8. %   smltn   simulation data structure
   
9. %   X   2D matrix i,j (states, time) NB this is the transpose of normal
   
10. %       return from the solver
    
11. %   S   3D matrix i,j,k (states, parameters, time)
    
12. 
    
13. % Validate the input/output arguments
    
14. if nargin~=2,
    
15.     error('Two input arguments are required');
    
16. end
    
17. if nargout<2 || nargout>3
    
18.     error('Two or three output arguments are required');
    
19. end
    
20. 
    
21. % Initialise the parameters
    
22. num_states = length(mdl.x_init);
    
23. num_param = length(mdl.theta);
    
24. 
    
25. % Reshape the initial parameter/sensitivity values for the ODE simulation)
    
26. xs_init = [mdl.x_init; reshape(mdl.S_init, num_states*num_param, 1)];
    
27. % Run the ODE simulation
    
28. [t, X_S] = feval(smltn.solver, mdl.ode, smltn.tspan, xs_init, [], mdl.theta);
    
29. 
    
30. num_time = length(t);
    
31. % Reshape the state matrix (states, time)
    
32. X = X_S(:,1:num_states)';
    
33. % Reshape the sensitivity derivatives (states, parameters, time)
    
34. if nargout==3,   
    
35.         S = reshape(X_S(:,num_states+1:end)', num_states, num_param, num_time);
    
36. end
    

复制代码

1. function dx = ODESenfun(t, x, k)
   
2. 
   
3. % ODESENFUN        Solve both system dynamic ODEs and Sensitivity Equations
   
4. % System Dynamic ODEs
   
5. %         x_dot = F*theta
   
6. %         S_dot = J*S + F
   
7. 
   
8. % Check input and output arguments are correctly supplied
   
9. if nargin~=3
   
10.     error('Three input arguments are required');
    
11. end
    
12. 
    
13. numStates = 10;
    
14. numPara = 11;
    
15. % Calculate the state derivatives
    
16. % Need to investigate a way to automatically generate these quantities (from a symbolic model)
    
17. 
    
18. %dx(1)=-k1*x(1)*x(6)+km1*x(2)+k4*x(4)-km4*x(1)*x(9)+k6*x(5);
    
19. %dx(2)=k1*x(1)*x(6) - km1*x(2) - k2*x(2) + km2*x(3)*x(7);
    
20. %dx(3)=k2*x(2) - km2*x(3)*x(7) - k3*x(3)*x(8) + km3*x(4) - k5W*x(3) + km5*x(5);
    
21. %dx(4)=k3*x(3)*x(8) - km3*x(4) - k4*x(4) + km4*x(1)*x(9);
    
22. %dx(5)=k5W*x(3) - km5*x(5) - k6*x(5);
    
23. %dx(6)=-k1*x(1)*x(6) + km1*x(2);
    
24. %dx(7)=k2*x(2) - km2*x(3)*x(7);
    
25. %dx(8)=-k3*x(3)*x(8) + km3*x(4);
    
26. %dx(9)=k4*x(4) - km4*x(1)*x(9);
    
27. %dx(10)=k6*x(5);
    
28. 
    
29. F = [-x(1)*x(6)      x(2)       0                   0              0           0          x(4)    -x(1)*x(9)         0        0            x(5);
    
30.      x(1)*x(6)      -x(2)      -x(2)             x(3)*x(7)         0           0          0          0               0        0             0;
    
31.         0             0         x(2)             -x(3)*x(7)    -x(3)*x(8)     x(4)        0          0             -x(3)     x(5)           0;
    
32.         0             0         0                   0           x(3)*x(8)     -x(4)      -x(4)     x(1)*x(9)         0        0             0;
    
33.         0             0         0                   0              0           0          0          0              x(3)     -x(5)         -x(5);
    
34.     -x(1)*x(6)       x(2)       0                   0              0           0          0          0               0        0             0;
    
35.         0             0        x(2)            -x(3)*x(7)          0           0          0          0               0        0             0;
    
36.         0             0         0                   0           -x(3)*x(8)    x(4)        0          0               0        0             0;
    
37.         0             0         0                   0              0           0         x(4)     -x(1)*x(9)         0        0             0;
    
38.         0             0         0                   0              0           0          0          0               0        0            x(5);];
    
39. 
    
40. dx = F*k;
    
41. 
    
42. % Calculate the sensitivity derivatives
    
43. % Need to investigate a way to automatically differentiate these quantities (from a symbolic model)
    
44. S = reshape(x(numStates+1:end), numStates, numPara);
    
45. 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;
    
46.       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;
    
47.          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;
    
48.       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;
    
49.          0                      0                k(9)                    0       -k(10)-k(11)      0        0            0           0          0;
    
50.     -k(1)*x(6)                 k(2)               0                      0          0         -k(1)*x(1)    0            0           0          0;
    
51.          0                     k(3)          -k(4)*x(7)                  0          0              0      -k(4)*x(3)     0           0          0;
    
52.          0                      0              -k(5)*x(8)                k(6)       0              0        0         -k(5)*x(3)     0          0;
    
53.     -k(8)*x(9)                  0                 0                      k(7)       0              0        0            0         -k(8)*x(1)   0;
    
54.          0                      0                 0                      0         k(11)           0        0            0           0          0;];
    
55.    
    
56. dS = J*S+F;
    
57. 
    
58. dx((numStates+1):(numStates*(numPara+1))) = reshape(dS,numStates*numPara,1);
    
59. 
    
60. 
    
61. 
    

复制代码

1. function S1 = S2DParam(S)
   
2. % S2DPARAM  Reshape the sensitivity matrix, S, to a 2D matrix
   
3. %   S - 3D matrix (states, param, time)
   
4. %   S1 - 2D matrix (states*time, param)
   
5. 
   
6. S1 = zeros(size(S,1)*size(S,3), size(S,2));
   
7. num_states = size(S,1);
   
8. num_time = size(S,3);
   
9. for i=1:num_states
   
10.     S1((i-1)*num_time+1:i*num_time, :) = shiftdim(S(i, :, :),1)';
    
11. end
    

复制代码

1. function smltn = smltnInit(solver, tspan)
   
2. % Initialize the parameters used in the ODE simulation
   
3. 
   
4. % Validate the input/output arguments
   
5. if nargin~=2
   
6.     error('2 input arguments must be supplied');
   
7. end
   
8. 
   
9. % Set the data structure's fields
   
10. smltn.solver = solver;
    
11. smltn.tspan = tspan;
    

复制代码

如果你有时间 不妨帮我看下， 在此我衷心感谢您的帮助与参与！！

ChaChing

1.水平有限, 2.时间有限
想想若是你, LZ会仔细看吗?
想法简化问题, 可能较容易...
看看 提问的智慧！！！！(发帖前请认真阅读)
http://forum.vibunion.com/forum/viewthread.php?tid=21991

lphtoby

回复 2 # ChaChing 的帖子

谢谢您的帮助 我会仔细学习的

ChaChing

1.总共6个程序, LZ想会有几个人帮忙看!?
2.LZ都知道plot(IKK, oed_a, r-)有问题了, 为何不检查矩阵变量的大小? ode_a为1*100, IKK为1*1, 当然画出直线!
3.个人水平专业有限, 仅能猜测修改, 在plot前加上IKK=IKK_normal*[1:100]/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.看看精华老帖, 一定有所得