Solve time-dependent ODE using the result from another time-dependent ODE

Hello everyone,

I needed to solve two ODEs of the following forms:

dRdS = 1 - B(S)*R(S) - A(S)*(R(S)^2); (1)

dWdS = - A(S)*R(S)*W(S) + R(S)*P(S); (2)

where B(S) , A(S) and P(S) are all deterministic functions depending on S, defined as follows:

Bs = linspace(78.809,80,100);

As = linspace(78.809,80,100);

Ps = linspace(78.809,80,100);

A = (2./(vm.*(As.^2))) .* ((sigma^2*vm)/(dv^2) + abs(alpha-beta*vm)/dv + r + 1/dtau);

B = -(2*(r-q)*Bs)./(vm*(Bs.^2));

P = (2./(vm.*Ps.^2)).*( ((-(sigma^2)*vm)/2)*(2*C/(dv^2)) - (abs(alpha-beta*vm)/2)*(2*C/dv) - (1/dtau)*C );

all the parameters above are defined beforehand.

I solved (1) by first writing a function dRdS:

function dRdS = dRdS(S,R,Bs,B,As,A)

B = interp1(Bs,B,S);

A = interp1(As,A,S);

dRdS = 1 - B.*R - A.*R*R;

end

and solved it using ode45:

S_span = [78.809 80];

IC = 0;

opts = odeset('RelTol',1e-2,'AbsTol',1e-4);

[S, R] = ode45(@(S,R) dRdS(S,R,Bs,B,As,A), S_span, IC, opts);

where I obtained a 53x2 matrix of [S, R].

I then moved on to write a function dWdS to solve (2):

function dWdS = dWdS(S,W,Bs,B,As,A,Ps,P,R)

B = interp1(Bs,B,S);

A = interp1(As,A,S);

P = interp1(Ps,P,S);

dWdS = -A.*R.*W-R.*P;

end

and using ode45:

IC_W = 0;

[S, W] = ode45(@(S,W) dWdS(S,W,Bs,B,As,A,Ps,P,Rm), S_span, IC_W);

While the ODE (1) was successfully solved, I encountered error message while solving ODE (2) that says:

"Error using odearguments (line 95)

@(S,W)DWDS(S,W,BS,B,AS,A,PS,P,R) returns a vector of length 53, but the length of initial

conditions vector is 1. The vector returned by @(S,W)DWDS(S,W,BS,B,AS,A,PS,P,R) and the

initial conditions vector must have the same number of elements."

In light of this error, I defined IC_W = zeros(53), but encountered another error: "Arrays have incompatible sizes for this operation."

So, I'd like to ask what exactly might have gone wrong in this implementation? If I'm adopting a wrong method for solving ODE (2), could someone enlighten me with a better method for solving such ODE which have the result from another ODE as part of the coefficient?

Sorry about the long question, and thanks so much!

[Torsten]

Why don't you solve both equations in one call to ODE45 ?

[Jan]

There are some abs() commands in the function to be integrated and interp1. Both functions produce non-smooth values. Matlab's ODE integrators are designed to handle smooth functions only. Running an integrator outside the conditions it is designed for, means running a random number generator with a very limited entropy. From the view point of numerical maths, this is no a scientific application of the algorithm. This does not cause the observed error.

The coorect way would be to let the integrator stop at each discontinuity using events, and to restart the integration.

We cannot run the posted code, because the variable vm is missing.

The output S and R of the first integration are not used for the second one.

Instead of using the pre-defined tables, it would be easier to move the definitions of A,B,P inside the function to be integrated. Torsten's suggestion to integrate both functions simultaneously, is a good idea also.

[Mingze Yin]

@Torsten Thanks for answering! Could you introduce to me a little more on how to solve both equations in one call? Sorry im actually not very familiar with the ode45 solver...

[Mingze Yin]

@Jan Thanks for answering! Allow me to spend some time understanding and implementing the ideas you shared :)

Answer 1

[Torsten]

function main
  IC = [0;0];
  S_span = [78.809 80];
  opts = odeset('RelTol',1e-4,'AbsTol',1e-4);
  [S, y] = ode45(@fun, S_span, IC, opts);
  R = y(:,1);
  W = y(:,2);
end

function dy = fun(S,y)
  R = y(1);
  W = y(2);
  vm = ...;
  sigma=...;
  dv = ...;
  alpha=...;
  beta=...;
  dtau=...;
  r=...;
  q=...;
  C=...;
  A = (2./(vm.*(S.^2))) .* ((sigma^2*vm)/(dv^2) + abs(alpha-beta*vm)/dv + r + 1/dtau);
  B = -(2*(r-q)*S)./(vm*(S.^2));
  P = (2./(vm.*S.^2)).*( ((-(sigma^2)*vm)/2)*(2*C/(dv^2)) - (abs(alpha-beta*vm)/2)*(2*C/dv) - (1/dtau)*C );
  dRdS = 1 - B.*R - A.*R*R;
  dWdS = -A.*R.*W-R.*P;
  dy = [dRdS;dWdS];
end

[Mingze Yin]

hey thanks so much for this, I tried this out and it worked fine! However, just one last follow-up question, your implementation seems to have the assumption that the s_span must be the same for both ODEs, but can this "solving multiple ODEs at the same time" method still work if the s_span is different for the ODEs to be solved?

[Torsten]

The ODE for W depends on R - so the s_span for W must be at least be a subset of the s_span of R.

So if the ODE for R is to be solved on [a1,b1] and the ODE for W is to be solved on [a2,b2], we have a1<=a2 and b1>=b2. Thus integrate the ODE for R alone on [a1,a2], then both ODEs together on [a2,b2] and the ODE for R again on [b2,b1].

[Mingze Yin]

Right that makes sense. Thanks for clarifying!

[Mingze Yin]

Hey sir sorry for bringing this question up again, but I encountered another issue while I was writing my code:

Building onto the previous ODEs that solve for R and W, now another ODE dVdS needed to be solved. I was trying to solve this dVdS together with the previous two since dVdS also depends on the values of R and W. However, there were some difficulties since I only know the terminal value (as opposed to the inital value) of dVdS, and so I couldnt figure out how to solve R, W and V altogether. I thought I could simply treat the terminal value as initial value but that would mean the s_span needs to be reversed? and reversing the s_span would not allow me to solve R and W I believe?

So I thought maybe I should solve dvds separately but then that would make R and W undefined. Is there a way to address this issue?

[Torsten]

Use bvp4c instead of ode45 or adjust the initial condition v_start for V at s=s_start such that you receive the desired terminal value v_terminal_desired for V at s=s_terminal. This can be done by using "fzero" to solve

v_terminal(v_start) - v_terminal_desired = 0.

[Pavitra Vaishali]

Hi I have same kind of problem, although I get martices as answer. Please, if anyone may help me with it