Solving coupled ode in MATLAB
显示 更早的评论
I am trying solve a system of coupled differential equantions (one of which is second order). I am searching for a numerical answer. However, code keeps giving me an error and I can't track it down. Would be greatful for any help! The code is below:
%function definitions
syms r(T) p(T);
ode1 = diff(r,T,2) == E^2 / ((m^2)*(c^2)) - (1- rs/r)*(c^2 + (h^2)/(r^2));
ode2 = diff(p,T) == (1/r^2) * L;
%function setup
f1= odeToVectorField(ode1);
F1= matlabFunction(f1);
f2= odeToVectorField(ode2);
F2= matlabFunction(f2);
odes = {F1,F2};
%initial conditions
r_0= 5*rs;
r_prime_0 = 5*rs;
r_double_prime_0 = 5*rs;
p_0 = 0;
p_prime_0 = PI/100;
conditions = [r_0;r_prime_0;r_double_prime_0;p_0;p_prime_0];
span = [0 100];
%solver
[rsol,psol] = ode45(odes,span,conditions);
采纳的回答
Star Strider
2020-1-27
This should get you started:
% %function definitions
syms r(T) p(T) c E h L m rs T Y
ode1 = diff(r,T,2) == E^2 / ((m^2)*(c^2)) - (1- rs/r)*(c^2 + (h^2)/(r^2));
ode2 = diff(p,T) == (1/r^2) * L;
% %function setup
[f1, Subs]= odeToVectorField([ode1 ode2]);
F1= matlabFunction(f1, 'Vars',{T,Y, c, E, h, L, m, rs})
% %initial conditions
r_0= 5*rs;
r_prime_0 = 5*rs;
r_double_prime_0 = 5*rs;
p_0 = 0;
PI = rand; % Provide Value
p_prime_0 = PI/100;
c = rand; % Provide Value
E = rand; % Provide Value
h = rand; % Provide Value
L = rand; % Provide Value
m = rand; % Provide Value
rs = rand; % Provide Value
conditions = [rand(3,1)]; % Three Equations = Three Initial Conditions
span = [0 100];
% %solver
[rsol,psol] = ode45(@(T,Y)F1(T,Y,c,E,h,L,m,rs),span,conditions);
figure
plot(rsol, psol)
grid
Note that if you want a numerical solution, the required parameter values must be defined first, before calling any of the numeric ODE solvers (such as ode45).
更多回答（1 个）
Pitt Lucht
2021-1-12
Hello,
I have a similar problem of three coupled differential equations. I have already adapted my code to fit the suggested code from Star Strider, but I'm still runnning into problems. The current error code I'm getting is "too many input arguments", my guess is there is a problem with the way I set up my initial conditions? If somebody has the time to have a quick look over my code it would be highly appreciated!
global kr kf lf dr dm p1 p2
% Parameter
kr = 1; % formation R-centers
kf = 0.1; % formation M-centers
lf = 100; % loss of F-centers
dr = 0.1; % decay R-center
dm = 1; % decay M-center
p1 = 10000;
p2 = 0;
% function definiton
syms R(T) M(T) F(T) T Y;
ode1 = diff(R,T) == -(dr*R) + (kr*M*F);
ode2 = diff(M,T) == (dr*R) - (dm*M) + (kf*(F^2)) - (kr*M*F);
ode3 = diff(F,T,2) == (dr*R) + (2*dm*M) - (kr*M*F) - (2*kf*(F^2)) - (lf*F)+p1;
% function setup
[f1, Subs] = odeToVectorField([ode1 ode2 ode3]);
F1 = matlabFunction(f1, 'Vars',{T, Y});
% initial states
R0 = 84.99;
M0 = 1.675;
F0 = 9.975;
conditions = [R0; M0; F0];
span = [0 10];
% solver
tic;
[Rsol, Msol, Fsol] = ode23s(@(T,Y)F1(T,Y,kr,kf,lf,dr,dm,p1),span,conditions);
toc;
%Log. plotting
loglog(Rsol, Msol, Fsol);
grid on;
5 个评论
Star Strider
2021-1-12
This should actually be a new Question.
That aside, this:
% initial states
R0 = 84.99;
M0 = 1.675;
F0 = 9.975;
conditions = [R0; M0; F0; 42]; % <— 4 DEs => 4 ICs! Added ‘42’ To Test Code
span = [0 10];
% solver
tic;
[Rsol, Msol] = ode23s(F1,span,conditions); % All Arguments Already Included In ‘F1’ After ‘subs’ Call
toc;
%Log. plotting
figure % Using 'log' Scale Not Appropriate With Negative Values, Use 'yyaxis' Instead
yyaxis left
plot(Rsol, Msol(:,[1 3 4])); % 2 Arguments Only!
yyaxis right
plot(Rsol, Msol(:,2))
grid on;
legend(string(Subs), 'Location','NW')
works.
Pitt Lucht
2021-1-12
Thank yo so much! And sorry for not opening a new question.
Star Strider
2021-1-12
My pleasure!
A Vote for my Answer would be appreciated!
Pitt Lucht
2021-1-12
It seems to not be possible to vote for comments, therefore I voted for your original answer to Lada Nuzhna since that helped out a lot too
Star Strider
2021-1-12
Thank you!
Yes, voting for an Answer is the only option in this instance.
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
