Solve time-dependent ODE using the result from another time-dependent ODE
2 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Mingze Yin
le 10 Déc 2021
Commenté : Pavitra Vaishali
le 12 Avr 2023
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 commentaires
Réponse acceptée
Torsten
le 12 Déc 2021
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 commentaires
Torsten
le 14 Déc 2021
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
le 12 Avr 2023
Hi I have same kind of problem, although I get martices as answer. Please, if anyone may help me with it
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!