MATLAB Answers

I am trying to solve 3 simultaneous ODE's. Which method to use?

3 views (last 30 days)
I am trying to solve 3 simultaneous ODE's. Which method should be used to solve these equations?? I have tried to solve it using ode45 method but I'm unable to get the answer.
dx1/df = F( f, x1, x2, x3 ) ]
dx2/df = F( f, x1, x2, x3 )
dx3/df = F( f, x1, x2, x3 )
f = [0.5 1]
x1(0.5)=0.185, x2(0.5)= 0.285, x3(0.5)= 0.53
gamma = 0.577;
xo1= 0.185; xo2= 0.285; xo3 =0.53;
theta=0.5;
q1=1; q2 = 0.317; q3 = 0.065;
F(x1,f) = ((q1)*(x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta)))-x1*(x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f)/(((x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f);
F(x2,f) = ((q2)*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))-x2*(x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f)/(((x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f);
F(x3,f)= ((q3)*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta)))-x3*(x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f)/(((x1-gamma*(x1*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(x2-gamma*(x2*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(x3-gamma*(x3*f-(xo3*(1-theta)))/(f-(1-theta))))*f);
  3 Comments
Dhiraj Deshmukh
Dhiraj Deshmukh on 6 Sep 2021
code i have used :-
function [dy] = countermulticompfcn(f,y)
theta = 0.5;
xo1 = 0.185;
xo2 = 0.285;
xo3 = 0.53;
gamma = 0.577;
q1=1;
q2 = 0.3543;
q3 = 0.3249;
dy = zeros(4,1)
dy(1) = ((q1)*(y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta)))-y(1)*(y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta)))))/(((y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta))))*f)
dy(2) = ((q2)*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))-y(2)*(y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta)))))/(((y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta))))*f)
dy(3) = ((q3)*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta)))-y(3)*(y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta)))))/(((y(1)-gamma*(y(1)*f-(xo1*(1-theta)))/(f-(1-theta)))+q2*(y(2)-gamma*(y(2)*f-(xo2*(1-theta)))/(f-(1-theta)))+q3*(y(3)-gamma*(y(3)*f-(xo3*(1-theta)))/(f-(1-theta))))*f)
end
solver:-
clear all
for index= 1:4
myprompt2=['Enter the initial value of Y:(',num2str(index),'):'];
initialval(index) = input(myprompt2);
end
lower = input('Enter Lower limit of Integration ');
upper = input('Enter Upper limit of Integration ');
[F,y] = ode45(@countermulticompfcn,[lower upper],initialval)
plot(F,y)
hold on

Sign in to comment.

Accepted Answer

Paul
Paul on 6 Sep 2021
This line:
dy = zeros(4,1);
Needs to be replaced with:
dy = zeros(3,1);
After that, note that the equation for dy(1) evaluates to NaN at f = 0.5 and to -Inf at f = 0 (I only explored between 0 and 1). In the original question it was stated that the the tspan input to ode45 is [0.5 1]. So that will cause a problem because the initial value of dy is NaN at f = 0.5 and that kills the entire solution. I ran from [0.51 1] and got a non-NaN result. In the follow up comment the user enters the tspan limits. If those limits cause the solver to hit either f = 0 or f = 0.5 exactly there will be a problem. So I guess the real question is if it's expected that the equations for dy have these singulariites and how they should be handled if the solver has to integrate through them.
  1 Comment
Dhiraj Deshmukh
Dhiraj Deshmukh on 7 Sep 2021
Thank you very much Paul. Problem is solved. Your suggestions really hepled.

Sign in to comment.

More Answers (1)

Steven Lord
Steven Lord on 6 Sep 2021
This documentation page includes two examples, the second of which you can use as a model for your code to solve a system of differential equations.

Community Treasure Hunt

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

Start Hunting!