Solving coupled ode in MATLAB

92 views (last 30 days)
Lada Nuzhna
Lada Nuzhna on 27 Jan 2020
Commented: Star Strider on 12 Jan 2021
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);

Accepted Answer

Star Strider
Star Strider on 27 Jan 2020
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).
  2 Comments
Lada Nuzhna
Lada Nuzhna on 27 Jan 2020
Thank you so so much!
Star Strider
Star Strider on 27 Jan 2020
As always, my pleasure!

Sign in to comment.

More Answers (1)

Pitt Lucht
Pitt Lucht on 12 Jan 2021
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 Comments
Pitt Lucht
Pitt Lucht on 12 Jan 2021
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
Star Strider on 12 Jan 2021
Thank you!
Yes, voting for an Answer is the only option in this instance.

Sign in to comment.

Categories

Find more on Symbolic Math Toolbox in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!