Может ли функция MATLAB ode45 принимать векторные массивы?

Question

Я получаю сообщение об ошибке в своем коде. Ошибка

" Error using odearguments (line 93) @(T,IP)V4.*(C+(1-ALPHA).*K4)./(C+K4)- 
(IR*IP) returns a vector of length 10000, but the length of initial 
conditions vector is 1. The vector returned by @(T,IP)V4.*(C+(1- 
ALPHA).*K4)./(C+K4)-(IR*IP) and the initial conditions vector must 
have the same number of elements " 

. Ошибка возникает из последней функции ODE45.

Это код:

clear all; clc; 
%-------------------------------------------------------------------------- 
% The simulation is based on the model described by DeYoung and Keizer in 
% the paper titled " A single-pool inositol 1,4,5-trisphosphate-receptor- 
% based model for agonist-stimulated oscillations in Ca2+ concentration" 
%-------------------------------------------------------------------------- 
%-------------------------------------------------------------------------- 
% Initial conditions 
%-------------------------------------------------------------------------- 
Ca_ER=10*10^-6; Ca_cyto=1.7*10^-6; 
Ir=1; alpha=.5; p_open3=0.15; 
%-------------------------------------------------------------------------- 
% Constants in micromolar 
%-------------------------------------------------------------------------- 
c0=4*10^-6; c1=.185; 
v1=6; v2=.11; v3=.09*10^-6; v4=1.2; 
k3=.1*10^-6; k4=1.1*10^-6; 
%-------------------------------------------------------------------------- 
% Receptor Binding Constants in micromolar per second 
%-------------------------------------------------------------------------- 
a1=400*10^-6; a2=0.2*10^-6; a3=400*10^-6; a4=0.2*10^-6; a5=20*10^-6; 
%-------------------------------------------------------------------------- 
% Receptor Dissociation Constants in micromolar 
%-------------------------------------------------------------------------- 
d1=0.13*10^-6; d2=1.049*10^-6; d3=.9434*10^-9; d4=.1445*10^-9; 
d5=.08234*10^-9; 
%-------------------------------------------------------------------------- 
% ODE describing Ca2+ concentrations in the cyctosol. Refer Ca2+ 
% oscillations 
%-------------------------------------------------------------------------- 
dc=@(t,c) (c1.*(v1.*(p_open3)+v2).*(Ca_ER)-c)-v3.*((c).^2)./(c.^2+(k3).^2); 
[t,c]=ode45(dc,linspace(0, 100, 10000),.15*10^-6); 
plot(t,c); 
%-------------------------------------------------------------------------- 
% Obtaining Ca_ER from the conservation condition. Refer Ca2+ oscillations 
%-------------------------------------------------------------------------- 
Ca_ER=(c0-c)./c1; 
figure(2); 
plot(t,Ca_ER); 
%-------------------------------------------------------------------------- 
%ODE describing IP3 production by Ca2+ feedback. Refer equation 4 
%-------------------------------------------------------------------------- 
dIP3= @(t,ip) v4.*(c+(1-alpha).*k4)./(c+k4)-(Ir*ip); 
[t,ip3]=ode45(dIP3,linspace(0, 100),.2*10^-6); 
plot(t,ip3);` 

Answer 1

Вы должны решить для обеих функций в то же время. Что-то вроде:

clear all; clc; 


%// Constants - units look nice as comments :-D 

Ca_ER = 10e-6; 
Ca_cyto = 1.7e-6; 
Ir  = 1; 
alpha = 0.5; 
p_open3 = 0.15; 

c0 = 4e-6; 
c1 = 0.185; 
v1 = 6; 
v2 = 0.11; 
v3 = 0.09e-6; 
v4 = 1.2; 
k3 = 0.1e-6; 
k4 = 1.1e-6; 

a1 = 400e-6; 
a2 = 0.2e-6; 
a3 = 400e-6; 
a4 = 0.2e-6; 
a5 = 20e-6; 

d1 = 0.13e-6; 
d2 = 1.049e-6; 
d3 = 0.9434e-9; 
d4 = 0.1445e-9; 
d5 = 0.08234e-9; 


%// Solving for f where: 
%// f(1) = Ca2+ conc 
%// f(2) = IP3 

df = @(t,f) [              ... 
    c1*(v1*p_open3 + v2)*Ca_ER - f(1) - v3*f(1)^2/(f(1)^2 + k3^2); ... 
    v4*(f(1)+k4*(1-alpha))/(f(1)+k4) - Ir*f(2)     ... 
]; 
f0 = [  ... 
    0.15e-6; ... 
    0.2e-6 ... 
]; 

[t,f] = ode45(df, linspace(0, 100, 10000), f0); 


%// Separate solutions 

c  = f(:,1); 
ip3 = f(:,2); 
Ca_ER1 =(c0-c)./c1; 


%// Plots 

figure(); plot(t, c); 
figure(); plot(t, Ca_ER1); 
figure(); plot(t, ip3);