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

1 次查看（过去 30 天）
Mingze Yin
Mingze Yin 2021-12-10
评论： Pavitra Vaishali 2023-4-12
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更早的评论
Mingze Yin
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
Mingze Yin 2021-12-12
@Jan Thanks for answering! Allow me to spend some time understanding and implementing the ideas you shared :)

请先登录，再进行评论。

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

采纳的回答

Torsten
Torsten 2021-12-12
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更早的评论
Torsten
Torsten 2021-12-14
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
Pavitra Vaishali 2023-4-12
Hi I have same kind of problem, although I get martices as answer. Please, if anyone may help me with it

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by