Solve time-dependent ODE using the result from another time-dependent ODE
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
Mingze Yin
am 10 Dez. 2021
Kommentiert: Pavitra Vaishali
am 12 Apr. 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 Kommentare
Akzeptierte Antwort
Torsten
am 12 Dez. 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 Kommentare
Torsten
am 14 Dez. 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
am 12 Apr. 2023
Hi I have same kind of problem, although I get martices as answer. Please, if anyone may help me with it
Weitere Antworten (0)
Siehe auch
Kategorien
Mehr zu Ordinary Differential Equations finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!