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

2 次查看（过去 30 天）

显示更早的评论

Mingze Yin 2021-12-10

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1607715-solve-time-dependent-ode-using-the-result-from-another-time-dependent-ode

评论： Pavitra Vaishali 2023-4-12

采纳的回答： Torsten

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!

4 个评论
显示 2更早的评论隐藏 2更早的评论

Mingze Yin 2021-12-12

@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 2021-12-12

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

请先登录，再进行评论。

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

采纳的回答

Torsten 2021-12-12

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1607715-solve-time-dependent-ode-using-the-result-from-another-time-dependent-ode#answer_853370

在 MATLAB Online 中打开

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

6 个评论
显示 4更早的评论隐藏 4更早的评论

Torsten 2021-12-14

在 MATLAB Online 中打开

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.