Numerical Solution of the System of four Coupled Nonlinear ODEs by Runge-Kutta Fourth Order Method

35 views (last 30 days)
I am trying to find out the Numerical Solution of the System of four Coupled Nonlinear ODEs using Runge-Kutta Fourth Order Method. I have developed the MatLab code for these equations. but i am not able to get correct solution. i have attached the Matlab code in this section. please any one notify me if any mistakes in this code..
a = 0; % Lower limits of integration
b = 1; % Upper limits of integration
h = 0.25; % h = delta_t
S(1) = 1;
E(1) = 0;
I(1) = 1;
R(1) = 0;
n = (b-a)/h ; % no of points
t = a + (0:n)*h
% Solution
func1 = @(S,I)-((0.5*I)+0.25)*S % function of integrate
func2 = @(E,I,S)(0.5*S*I)-E % function of integrate
func3 = @(I,E) (0.75*E)-(I) % function of integrate
func4 = @(R,I) (0.75*I)-(0.25*R) % function of integrate
% func_exact = @(t) ;
for i = 1:n
k1 = func1( I(i),S(i));
k2 = feval(func1, I(i)+(k1/2),S(i)+(k1/2));
k3 = feval(func1, I(i)+(k2*h/2),S(i)+(k2*h/2));
k4 = feval(func1, I(i)+(k3*h),S(i)+(k3*h));
l1 = feval(func2, S(i),I(i),E(i));
l2 = feval(func2, S(i)+(l1*h/2),I(i)+(l1*h/2),E(i)+(l1*h/2));
l3 = feval(func2, S(i)+(l2*h/2),I(i)+(l2*h/2),E(i)+(l2*h/2));
l4 = feval(func2, S(i)+(l3*h),I(i)+(l3*h),E(i)+(l3*h));
m1 = feval(func3, E(i),I(i));
m2 = feval(func3, E(i)+(m1*h/2),I(i)+(m1*h/2));
m3 = feval(func3, E(i)+(m2*h/2),I(i)+(m2*h/2));
m4 = feval(func3, E(i)+(h*m3),I(i)+(m3*h));
n1 = feval(func4, I(i),R(i));
n2 = feval(func4, I(i)+(n1*h/2),R(i)+(n1*h/2));
n3 = feval(func4, I(i)+(n2*h/2),R(i)+(n2*h/2));
n4 = feval(func4, I(i)+(n3*h),R(i)+(n3*h));
S(i+1) = S(i) + (h/6)*(k1+(2*(k2+k3))+k4)
E(i+1) = E(i) + (h/6)*(l1+(2*(l2+l3))+l4)
I(i+1) = I(i) + (h/6)*(m1+(2*(m2+m3))+m4)
R(i+1) = R(i) + (h/6)*(n1+(2*(n2+n3))+n4)
end
format bank
disp('Numerical results'),disp(' t S E I R');
disp ([t',S',E',I', R'])

Answers (1)

James Tursa
James Tursa on 30 Jul 2020
Edited: James Tursa on 30 Jul 2020
I will point out some problems, and then suggest a much easier way to do this.
Start with your first couple of lines in the loop:
k1 = func1( I(i),S(i));
k2 = feval(func1, I(i)+(k1/2),S(i)+(k1/2));
func1 is dS/dt, so that S(i)+(k1/2) part is correct (almost). But the I(i)+(k1/2) part is not correct because k1 is dS/dt and it should be dI/dt for this part. And you are missing the h. That is, these lines should look something like this instead:
k1 = func1( I(i),S(i)); % dS/dt at t
m1 = feval(func3, E(i),I(i)); % dI/dt at t
k2 = feval(func1, I(i)+(m1/2)*h,S(i)+(k1/2)*h); % dS/dt at t+h/2
You have similar problems with the rest of your code. You need to evaluate all of the k1, l1, m1, n1 terms first. Then use those to evaluate all of the k2, l2, m2, n2 terms. Etc. Don't mix the k, l, m, n variables in your calculations. k's go with S's, l's go with E's, m's go with I's, and n's go with R's. And you need to include h in these calculations.
Having said all of that, this would be much easier and less code to write if you vectorized everything and then wrote vector equations for your derivative functions and RK4 code. That way you would not have to duplicate all the RK4 stuff. E.g., define a y vector as
y(1) = S
y(2) = E
y(3) = I
y(4) = R
Then write one derivative function that takes a 4-element y as an input and returns a 4-element dy/dt as an output. That way you only need one set of RK4 code instead of four sets. Another benefit of coding it this way is you can pass the function handle and initial conditions directly to the ode45( ) function to double check your answer.

Categories

Find more on Dates and Time 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!