MATLAB Answers

I dont know what to do

1 view (last 30 days)
Meor Hasan Meor Jumat
Meor Hasan Meor Jumat on 12 Jun 2021
Commented: Walter Roberson on 16 Jun 2021
% Fourth Order Runge Kutta Method to solve system of ODEs
% y1'= f1(t, y1, y2, y3, y4);
% y2'= f2(t, y1, y2, y3, y4);
% y3'= f2(t, y1, y2, y3, y4);
% y4'= f2(t, y1, y2, y3, y4);
% Fourth Order RK method for solving ODE
% y' = f(t, y)
% y_n+1 =y_n + (h/6) * (delta_y1 + 2*delta_y2 + 2*delta_y3 + delta_y4)
% delta_y1 = f(tn,yn)
% delta_y2 = f(tn + h/2, yn + h*(delta_y1/2))
% delta_y3 = f(tn + h/2, yn + h*(delta_y2/2))
% delta_y4 = f(tn + h, yn + h*delta_y3)
% dy1/dx = -*y1*(y2+y3)
% dy2/dx = 0.47*(1.41*10^-9)*y1*(y2+y3)-0.1*y2
% dy3/dx = 0.53*(1.41*10^-9)*y1*(y2+y3)-0.1*y3
% dy4/dx = 4 - 0.3 * y2 - 0.1 * y1
% Numerically integrate above from x = 0 to x = 60 using step size of 0.001
% IC: y1(x=0) = 8.18*10**7; y2(x=0) = 5.0604*10**3; y3(x=0) = 4.4876*10**3; y4(x=0) = 4590;
% y1exact = ;
% y2exact = ;
% y3exact = ;
% y4exact = ;
clear
clc
clf
% Initial Value and time step
a = 1.41*10^-9; % Lower limit of Integration
b = 2*10^-9; % Upper limit of Integration
h = 0.001; % h = delta_x
y1(1) = 8.18*10^7; % S(1)
y2(1) = 5.0604*10^3; % I_1(1)
y3(1) = 4.4876*10^3; % I_2(2)
y4(1) = 4590; % R(1)
n = (b-a)/h + 1; % no. of points
% Solution
x(1) = 0;
func1 = @(y1, y2, y3, y4) (-1.41*10^-9 * y1 * (y2 + y3)); % function to integrate dS(t)/dx
func2 = @(y1, y2, y3, y4) ((0.45 * 1.41*10^-9 * y1 * (y2 + y3)) - ((1/10) * y2)); % function to integrate dI_1(t)/dx
func3 = @(y1, y2, y3, y4) (((1 - 0.45) * 1.41*10^-9 * y1 * (y2 + y3)) - (((1/15) + (1/50)) * y3)); % function to integrate dI_2(t)/dx
func4 = @(y1, y2, y3, y4) (((1/10) * y2) + ((1/15) * y3)); % function to integrate dR(t)/dt
% func_exact = @(x) ;
for i = 1:n-1
x(i+1) = 1*h;
delta_y1_1 = h * feval(func1, y1(i));
delta_y2_1 = h * feval(func1, y1(i) + (delta_y1_1)/2);
delta_y3_1 = h * feval(func1, y1(i) + (delta_y2_1)/2);
delta_y4_1 = h + feval(func1, y1(i) + delta_y3_1);
y1(i+1) = y1(i) + (1/6) + (delta_y1_1 + 2*delta_y2_1 + 2*delta_y3_1 + delta_y4_1);
delta_y1_2 = h * feval(func2, y1(i), y2(i));
delta_y2_2 = h * feval(func2, y1(i) + (delta_y1_1)/2, y2(i) + (delta_y1_2)/2);
delta_y3_2 = h * feval(func2, y1(i) + (delta_y2_1)/2, y2(i) + (delta_y2_2)/2);
delta_y4_2 = h + feval(func2, y1(i) + delta_y3_1, y2(i) + delta_y3_2);
y2(i+1) = y2(i) + (1/6) + (delta_y1_2 + 2*delta_y2_2 + 2*delta_y3_2 + delta_y4_2);
delta_y1_3 = h * feval(func3, y1(i), y2(i), y3(i));
delta_y2_3 = h * feval(func3, y1(i) + (delta_y1_1)/2, y2(i) + (delta_y1_2)/2, y3(i) + (delta_y1_3)/2);
delta_y3_3 = h * feval(func3, y1(i) + (delta_y2_1)/2, y2(i) + (delta_y2_2)/2, y3(i) + (delta_y2_3)/2);
delta_y4_3 = h + feval(func3, y1(i) + delta_y3_1, y2(i) + delta_y3_2, y3(i) + delta_y3_3);
y3(i+1) = y3(i) + (1/6) + (delta_y1_3 + 2*delta_y2_3 + 2*delta_y3_3 + delta_y4_3);
delta_y1_4 = h * feval(func4, y1(i), y2(i), y3(i), y4(i));
delta_y2_4 = h * feval(func4, y1(i) + (delta_y1_1)/2, y2(i) + (delta_y1_2)/2, y3(i) + (delta_y1_3)/2, y4(i) + (delta_y1_4)/2);
delta_y3_4 = h * feval(func4, y1(i) + (delta_y2_1)/2, y2(i) + (delta_y2_2)/2, y3(i) + (delta_y2_3)/2, y4(i) + (delta_y2_4)/2);
delta_y4_4 = h + feval(func4, y1(i) + delta_y3_1, y2(i) + delta_y3_2, y3(i) + delta_y3_3, y4(i) + delta_y3_4);
y4(i+1) = y4(i) + (1/6) + (delta_y1_4 + 2*delta_y2_4 + 2*delta_y3_4 + delta_y4_4);
end
x;
y1;
y2;
y3;
y4;
  3 Comments
Meor Hasan Meor Jumat
Meor Hasan Meor Jumat on 12 Jun 2021
the results doesnt come out

Sign in to comment.

Answers (1)

Jan
Jan on 12 Jun 2021
Edited: Jan on 12 Jun 2021
Remove the semicolons from these lines:
x;
y1;
y2;
y3;
y4;
The trailing semicolon suppresses the output. Therefore "x;" shows x without showing x. Of course you do not see anything. You will see, that all variables are scalars, because you calculate 0 steps
a = 1.41*10^-9; % Lower limit of Integration
b = 2*10^-9; % Upper limit of Integration
h = 0.001; % h = delta_x
n = (b - a) / h + 1
Now n = 1.00000059
Then:
for i = 1:n-1
performs no iterations and the trajectory stays at the initial values.
By the way: -1.41*10^9 is an expensive power operations, while -1.41e9 is a cheap constant.
I'd prefer to convert:
func4 = @(y1, y2, y3, y4) (((1/10) * y2) + ((1/15) * y3));
% to
func4 = @(y1, y2, y3, y4) y2 / 10 + y3 / 15;
There is something wrong with your Runge-Kutta code. Why does the 2nd component depend on the result of the 1st one, but not the other way around?
func1 = @(y1, y2, y3, y4), so it depends on 4 inputs. But here you provide only 1 input:
delta_y1_1 = h * feval(func1, y1(i));
This bug did not cause problems yet, because the loop is not entered at all.
Calculating the functions elementwise demands for writing a Runge Kutta code for each number of dimensions. It is much easier to use vectors instead:
y = zeros(n, 4);
x = zeros(n, 1);
y(1, :) = [8.18e7, 5.0604e3; 4.4876e3; 4590];
x(1) = a; % Not 0
fcn = @(x, y) [-1.41e-9 * y(1) * (y(2) + y(3)), ...
0.45 * 1.41e-9 * y(1) * (y(2) + y(3)) - y(2) / 10, ...
0.55 * 1.41e-9 * y(1) * (y(2) + y(3)) - y(3) * 13 / 150, ...
y(2) / 10 + y(3) / 15];
for i = 1:n-1
x(i + 1) = x(i) + h; % Not: "x(i+1) = 1*h", but maybe: "i * h"
d1 = h * fcn(x(i), y(i, :));
d2 = h * fcn(x(i), y(i, :) + d1 / 2);
d3 = h * fcn(x(i), y(i, :) + d2 / 2);
d4 = h + fcn(x(i), y(i, :) + d2);
y1(i+1, :) = y1(i, :) + (d1 + d2 + d2 + d4) / 6;
end
Do you see the benefit of using vectors? Less clutter, less bugs.
Although your fcn does not depened on x, I've inserted it here such that you can use the Runge Kutta part in general.
  8 Comments
Walter Roberson
Walter Roberson on 16 Jun 2021
Using the hint from James, let us experiment:
% Initial Value and time step
a = 21; % Lower limit of Integration
b = 31; % Upper limit of Integration
h = 0.5; % h = delta_x
n = (b-a)/h + 1; % no. of points
n = 5;
y = zeros(n, 4, 'sym');
x = zeros(n, 1, 'sym');
y(1, :) = [10276; 751; 11; 9395];
x(1) = a; % Not 0
fcn = @(x, y) [-10219 * y(1) * (y(2) + y(3)),...
1.26 * 10219 * y(1) * (y(2) + y(3)) - 7.31 * y(2), ...
-0.26 * 10219 * y(1) * (y(2) + y(3)) - (0.11 + 1.26) * y(3), ...
7.31 * y(2) + 0.11 * y(3) ]; % system function
for i = 1:n-1
x( i + 1 ) = x(i) + h; % Not: "x(i+1) = 1*h", but maybe: "i * h"
d1 = h * fcn(x(i), y(i, :));
d2 = h * fcn(x(i), y(i, :) + d1 / 2);
d3 = h * fcn(x(i), y(i, :) + d2 / 2);
d4 = h + fcn(x(i), y(i, :) + d2);
y(i+1, :) = y(i, :) + (d1 + 2*d2 + 2*d3 + d4) / 6;
end
x
x = 
y
y = 
vpa(y)
ans = 
plot ( x, y )
This suggests that there are other problems with the system.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!