# I dont know what to do

1 ビュー (過去 30 日間)
Meor Hasan Meor Jumat 2021 年 6 月 12 日
コメント済み: Walter Roberson 2021 年 6 月 16 日
% 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 件のコメント表示非表示 2 件の古いコメント
Meor Hasan Meor Jumat 2021 年 6 月 12 日
the results doesnt come out

サインインしてコメントする。

### 回答 (1 件)

Jan 2021 年 6 月 12 日

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 件のコメント表示非表示 7 件の古いコメント
Walter Roberson 2021 年 6 月 16 日
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.

サインインしてコメントする。

### Community Treasure Hunt

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

Start Hunting!

Translated by