MATLAB Answers

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

조회 수: 3(최근 30일)
Dhiraj Deshmukh
Dhiraj Deshmukh 2021년 9월 6일
댓글: Dhiraj Deshmukh 2021년 9월 7일
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
Dhiraj Deshmukh
Dhiraj Deshmukh 2021년 9월 6일
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

댓글을 달려면 로그인하십시오.

채택된 답변

Paul
Paul 2021년 9월 6일
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
Dhiraj Deshmukh
Dhiraj Deshmukh 2021년 9월 7일
Thank you very much Paul. Problem is solved. Your suggestions really hepled.

댓글을 달려면 로그인하십시오.

추가 답변(1개)

Steven Lord
Steven Lord 2021년 9월 6일
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!

Translated by