How to call multiple functions in ode 45?

Question

gorilla3 2017-11-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/367015-how-to-call-multiple-functions-in-ode-45

评论： gorilla3 2017-11-23

I'm trying to solve 2 systems of differential equations using ode45, however I get the error code:

Error using odearguments (line 21) When the first argument to ode45 is a function handle, the tspan argument must have at least two elements.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);
Error in CBF_2006_f_v1 (line 49)
    [ABP,P1,P2,xq,xm,xc,xm1,Cav] = ode45(@vsa, @aut, tspan, var0);

There are 2 systems: 4 equations each. The unknowns are (V_sa,P1,P2,xq,xm,xc,xm1,Cav), all the rest are parameters that have been declared at the beginning of the code.

The code main code that recalls the functions is:

%Parameter declaration .... (Rsa,V_sa_b,P_ic ...)
%%Diff eq system
      ABP= linspace(40,170,131);
      delta=deltaCa_p;
      for i=1:1:length(ABP)
      tspan=[0 100];
      var0=[12; 97.6; 49.67; 0; 0; 0; 0; 0.205];
      [ABP,P1,P2,xq,xm,xc,xm1,Cav] = ode45(@vsa, @aut, tspan, var0);
      end

The 2 functions are stored in separate files and their code is:

function dvdt = vsa(ABP,var,delta,R_la,R_sa_b,V_sa_b,R_sv,P_ic,P_v,k_ven,P_v1)
dvdt = [(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1)))*(var(2)-P_ic) + var(8)*(( (ABP(i) -var(2))/(R_la+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-var(3))/(R_sv+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-P_ic)*(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))))  /var(8)) ;
    ( (ABP(i) -var(2))/(R_la+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-var(3))/(R_sv+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-P_ic)*(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))))  /var(8);
                  ((var(2)-var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) +R_sv)-(var(3)-P_v)/R_lv)*1/(1/(k_ven*(var(3)-P_ic-P_v1))) ;
                      0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))] ;
V_sa= var(:,1);
P1= var(:,2);
P2 = var(:,3);
Cav=var(:,8);
xq= var(:,4);
xc= var(:,5);
xm= var(:,6);
xm1= var(:,7);
dvdt=[V_sa;P1;P2;Cav; xq;xc;xm;xm1;Cav];
end

and

    function dxdt= aut(ABP,var,G_q,R_sa_b,V_sa_b,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)
    dxdt=[                                     (-var(4)+G_q*( ( (var(2)- var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) + R_sv) ) -q_b /q_b) )/tau_q ;
                                               (-var(5) +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
                                                                  var(7);
                                                                     (eps*u-tau2*var(7)-var(6))/tau1^2 ];
    V_sa= var(1);
    P1= var(2);
    P2 = var(3);
    xq= var(4);
    xc= var(5);
    xm= var(6);
    xm1= var(7);
    dxdt=[P1;P2;xq;xc;xm;xm1];
    end

Could you please help me find the mistakes?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2017-11-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/367015-how-to-call-multiple-functions-in-ode-45#answer_291104

在 MATLAB Online 中打开

I would nest them into one function:

bothfcns = @(ABP, var) [@(ABP,var) vsa(ABP,var,delta,R_la,R_sa_b,V_sa_b,R_sv,P_ic,P_v,k_ven,P_v1); @(ABP,var) aut(ABP,var,G_q,R_sa_b,V_sa_b,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)];

Then use ‘bothfcns’ as the function argument to ode45.

NOTE — I cannot run your code to test this, so I am listing this as UNTESTED CODE.

29 个评论
显示 27更早的评论隐藏 27更早的评论

Star Strider 2017-11-14

在 MATLAB Online 中打开

My pleasure.

I created ‘bothfcns’ to vertically concatenate ‘vsa’ and ‘aut’ so they would work in ode45. It is only a function handle, and it is not possible to get anything else from it, especially the result of the ode45 integration. Those results are in the independent variable ‘ABP’ and the dependent variable matrix ‘var’.

The way I read your code:

xq=var(:,12);
xc=var(:,13);
xm=var(:,14);

(you need to verify that the column references are correct) are in the integrated output after ode45 completes.

I have absolutely no idea what ‘xm_res’ and the rest are, since they appear in your last comment for the first time. However to sum the corresponding columns in ‘var’, just select the columns and sum across the rows:

sum_xqxcxm = sum(var(:, 12:14), 2);

(again assuming I have the columns correct). The result, ‘sum_xqxcxm’ will be a column vector. The last element is simply ‘sum_xqxcxm(end)’. Alternatively, if you only want to sum the last element of each vector, that would be:

sum_xqxcxm = sum(var(end, 12:14), 2);

It is neither appropriate nor likely possible to do the column summation while ode45 is calculating them. You have to wait until it is finished, and then do the calculations on the results.

gorilla3 2017-11-15

在 MATLAB Online 中打开

Hi,

this is the code I used:

sum_xqxcxm = sum(var(end, 12:14), 2);

Why would var be a set of random numbers?it should be the vector of solutions of the ode45...

Also, as you suggested, I inserted 'bothfncs' as the function argument to ode45 and this is the error message:

Error using feval
Undefined function or variable 'bothfncs'.

For completion here's the full code:

%%Parameters
R_la= 0.4; % mmHg s/ml
R_sa_b= 5.03;
R_sv= 1.32;
R_lv= 0.56;
P_v= 6;
V_la=1;
V_sa_b= 12;
P_ic= 10;
k_ven= 0.186; %[/ml] elastance coeff for venous compliance
P_v1= -2.25; % pressure offset for venous compl
V_vn= 28; %offset of vG_q= 3; %%gain of flow based feedback mechanism
tau_q= 20; %[s] time const of flow based feedback mechanism
Pa_co2_b= 40; % baseline arterial CO2 pressure
tau_co2= 40; %time const for arterial co2 p
tau1= 2; %1st time const of neural feedback
tau2= 4;
tau_g= 1; % t const for tissue O2 concentration
C_a_p=2.87; %ampl of + change in art comp [ml/mmHg]
C_a_n= 0.164;
g=0.2; %baseline non dim tissue O2 concentration
E=0.4; %baseline O2 extraction fraction
K= 0.15; %scaling between activation and O2 demand
V0= 0.02; %scaling factor for BOLD response
q_b=12.5;
G_q=3; %gain of flow based feedback mechanism [ml/mmHg]
Pa_co2=40; 
Ca_b= 0.205;
eps=1;
u=0.5; %0 for steady state
Pa_b=100;
ka=3.68;
deltaCa_p=2.87;
deltaCa_n=0.164;
%%Diff eq system
     ABP= linspace(40,170,131);  
      delta=deltaCa_p;
      tspan=[0 100];
      var0=[12; 97.6; 49.67; 0; 0; 0; 0; 0.205];
      bothfcns = @(ABP, var) [@(ABP,var) vsa(ABP,var,delta,R_la,R_sa_b,V_sa_b,R_sv,P_ic,P_v,k_ven,P_v1); @(ABP,var) aut(ABP,var,G_q,R_sa_b,V_sa_b,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)];
  [ABP,var] = ode45('bothfncs',tspan, var0);
sum_xqxcxm = sum(var(:, 12:14), 2);
function dvdt = vsa(ABP,var,delta,R_la,R_sa_b,V_sa_b,R_sv,P_ic,P_v,k_ven,P_v1)
      for i=1:1:length(ABP)
  dvdt = [(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1)))*(var(2)-P_ic) + var(8)*(( (ABP(i) -var(2))/(R_la+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-var(3))/(R_sv+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-P_ic)*(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))))  /var(8)) ;
      ( (ABP(i) -var(2))/(R_la+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-var(3))/(R_sv+0.5*(R_sa_b*(V_sa_b/var(1))^2)) - (var(2)-P_ic)*(0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))))  /var(8);
                    ((var(2)-var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) +R_sv)-(var(3)-P_v)/R_lv)*1/(1/(k_ven*(var(3)-P_ic-P_v1))) ;
                        0.5*delta*(1- ((cosh(4*(var(6)+var(5)-var(4))/delta))-1)/(cosh(4*(var(6)+var(5)-var(4))) +1))] ;
V_sa= var(:,1);
P1= var(:,2);
P2 = var(:,3);
Cav=var(:,8);
xq= var(:,4);
xc= var(:,5);
xm= var(:,6);
xm1= var(:,7);
dvdt=[V_sa;P1;P2;Cav; xq;xc;xm;xm1;Cav];
end
end
function dxdt= aut(ABP,V_sa,P1,P2,G_q,R_sa_b,V_sa_b,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)
dxdt=[                                     (-var(4)+G_q*( ( (var(2)- var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) + R_sv) ) -q_b /q_b) )/tau_q ;
                                           (-var(5) +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
                                                              var(7);
                                                                 (eps*u-tau2*var(7)-var(6))/tau1^2 ];
V_sa= var(1);
P1= var(2);
P2 = var(3);
xq= var(4);
xc= var(5);
xm= var(6);
xm1= var(7);
dxdt=[P1;P2;xq;xc;xm;xm1];
end

Star Strider 2017-11-15

在 MATLAB Online 中打开

I managed to get ‘vsa’ working by using single subscripts only to refer to elements of ‘var’.

Note that both ‘dvdt’ and ‘dxdt’ return (4x1) column vectors.

In ‘aut’, you use several variables you apparently have not defined, so even when I add them to the argument list, the function fails with a ‘Not enough input arguments’ error. You also did not include ‘var’ as the second argument to it, as required, Adding that throws other errors when I tried to run it independently.

As an initial approach, get each of your two differential equations working with ode45 first. When they run, try them with my original ‘bothfcns’ anonymous function.

Your code is so complicated that I can’t even begin to figure it out. I will help you with it as I can. You have to find out what problems you are having with it.

Also, my original function vertical concatenation idea works correctly in this test code:

f1 = @(x) [x(1).^2; x(2).^3];
f2 = @(x) [sqrt(x(1)); x(1)*x(2)];
f3 = @(x) [f1(x); f2(x)];
x = [2; 3];
y = f3(x);

I have no idea what the problem was when you were unable to use my ‘bothfcns’ anonymous function. I cannot test it with them because I cannot get ‘aut’ to work.

gorilla3 2017-11-15

编辑：gorilla3 2017-11-15

在 MATLAB Online 中打开

I would like to solve the 2 systems of equations, then (from aut) extrapolate the solutions for xc,xm,xq and sum them. Subsequently I would like to check whether the sum is positive or negative and according to the sign, allocate a value to the parameter delta. If the sum is + then delta=2.8, else delta=0.16. Delta is used in system vsa.

I fixed the function aut, and re-named it "second"(because I read file name and function name should not be the same), so it's working now:

function dxdt= second(ABP,var, V_sa,P1,P2,Cav,xq,xm,xc,xm1,R_sa_b,V_sa_b,R_sv,Pa_co2,Pa_co2_b,tau_co2,eps,u,tau2,tau1)
dxdt=[                                     (-var(4)+G_q*( ( (var(2)- var(3))/(0.5 *(R_sa_b*(V_sa_b/var(1))^2) + R_sv) ) -q_b /q_b) )/tau_q ;
                                           (-var(5) +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
                                                              var(7);
                                                                 (eps*u-tau2*var(7)-var(6))/tau1^2 ];
V_sa= var(1);
P1= var(2);
P2 = var(3);
xq= var(4);
xc= var(5);
xm= var(6);
xm1= var(7);
dxdt=[V_sa;P1;P2;xq;xc;xm;xm1];
end

__________________

Thank you Greg, but that's not the issue.

gorilla3 2017-11-22

在 MATLAB Online 中打开

Hi, I improved and simplified the code. Here it is:

clear all 
clc
%%Diff eq system
function second
[ABP,par] = ode45('first',[0 100], [12 97.6 49.67 0.205 0 0 0 0]);
sum_xqxcxm = sum(par(:, 4:6), 2);
function dvdt = first(ABP,par)
    %%Parameters
R_la= 0.4; % mmHg s/ml
R_sa_b= 5.03;
R_sv= 1.32;
R_lv= 0.56;
P_v= 6;
V_la=1;
V_sa_b= 12;
P_ic= 10;
k_ven= 0.186;
P_v1= -2.25;
V_vn= 28; 
tau_q= 20; 
Pa_co2_b= 40; 
tau_co2= 40; 
tau1= 2; 
tau2= 4;
tau_g= 1; 
C_a_p=2.87; 
C_a_n= 0.164;
g=0.2; 
E=0.4; 
K= 0.15; 
V0= 0.02; 
q_b=12.5;
G_q=3; 
Pa_co2=40; 
Ca_b= 0.205;
eps=1;
u=0.5; 
Pa_b=100;
ka=3.68;
deltaCa_p=2.87;
deltaCa_n=0.164;
 ABP= linspace(40,170,131); 
       for i=1:1:length(ABP)
% %calculate XSUM first and find 'i' values
% if i== (1||2||n...)
% delta=deltaCa_p; 
% elseif i==(4||5...100)
% delta=deltaCa_n;
% end
delta=2.15 % this is just to make it run but needs to be fixed, see explanation below
dV_sa= dCa*(par(2)-P_ic) + par(8)*dP1;
dP1= 1/par(8) * ((ABP(i)-par(2))/(R_la+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) - (par(2)-par(3))/(R_sv+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) -dCa*(par(2)-P_ic));
dP2=1/(1/(k_ven*(par(3)-P_ic-P_v1)))*((par(2)-par(3))/(R_sv+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) -(par(3)-P_v)/R_lv);
dCa=0.5*delta*(1- ((cosh(4*(par(6)+par(5)-par(4))/delta))-1)/(cosh(4*(par(6)+par(5)-par(4))) +1));
dxq= (-par(4)+G_q*( ( (par(2)- par(3))/(R_sv+0.5*(R_sa_b/(par(1)/V_sa_b)^2)) ) -q_b /q_b) )/tau_q ;
dxc=(-par(5) +0.3+3*tanh(Pa_co2/Pa_co2_b -1.1))/tau_co2;
dxm1=par(7);
dxm=(eps*u-tau2*par(7)-par(6))/tau1^2;
V_sa= par(1);
P1= par(2);
P2 = par(3);
Cav=par(8);
xq= par(4);
xc= par(5);
xm= par(6);
xm1= par(7);
par=[V_sa;P1;P2;Cav;xq;xc;xm;xm1];
dvdt=[dV_sa;dP1;dP2;dCa;dxq;dxc;dxm1;dxm];
end
end
end

As you can see, at each time step the value of ABP changes (for i=1:1:length(ABP)). I would now like to extract the solutions for xq,xc,xm at each time step t. Then I would like to sum these 3 solutions for each time step and print them in an array.

sum_xqxcxm = sum(par(:, 4:6), 2);

Then check for each value int he array if the value is >0 or <0. Assign delta=2.87 if positive and delta=0.16 if negative. Where delta is a parameter in the system of equations.

Torsten 2017-11-23

在 MATLAB Online 中打开

If you save the complete code in one file, the first part also has to be included in a function. So start your file with

function main
ABP= linspace(40,170,131);  
for i=1:length(ABP)
  [T,Y] = ode45(@(t,y)ABP(t,y,ABP(i)),[0 100], [12 97.6 49.67 0.205 0 0 0 0]);
end
function dvdt = first(t,par,abp)
...

and name your file "main.m".

Best wishes

Torsten.

gorilla3 2017-11-23

Doing that gives me the error: "unable to define function 'main' because it has the same name as the file"

请先登录，再进行评论。

How to call multiple functions in ode 45?

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

29 个评论
显示 27更早的评论隐藏 27更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

How to call multiple functions in ode 45?

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

29 个评论 显示 27更早的评论隐藏 27更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

29 个评论
显示 27更早的评论隐藏 27更早的评论