RK 4 method for ODES

1 次查看(过去 30 天)
Syed Abdul Rafay
Syed Abdul Rafay 2022-12-15
When I run this code it gives me error. I tried to change some value but no luck
function RK4(a,b,N)
if nargin<3, a=0; b=1; N=20; end
h=(b-a)/N;
y(1)=1;
y1(1)=2;
for i=1:N
t=(i-1)*h;
k1=f(t,y(i),y1(i));
k2=f(t+h/2,y(i)+k1*h/2,y1(i)+k1*h/2);
k3=f(t+h/2,y(i)+k2*h/2,y1(i)+k2*h/2);
k4=f(t+h,y(i)+k3*h,y1(i)+k3*h);
y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4)*h;
end
for i=1:N+1
t=(i-1)*h;
exact(i)=exp(t)+t;
end
error=max(abs(exact-y))
disp(y(N));
figure;
plot(i,y,'LineWidth',2);
xlabel('X');
ylabel('Y');
end
function fty=f(t,y,y1)
fty=-y1+exp(-7*t)*(t-y)+2*exp(t)+exp(-6*t)+1;
end
Output
Index exceeds the number of array elements. Index must not exceed 1.
Error in RK4 (line 28)
k1=f(t,y(i),y1(i));

回答(1 个)

Steven Lord
Steven Lord 2022-12-15
y(1)=1;
y1(1)=2;
for i=1:N
At the start of the loop (during iteration 1) y and y1 each have 1 element.
t=(i-1)*h;
k1=f(t,y(i),y1(i));
At iteration 1 this line uses element 1 of the y and y1 variables.
At iteration 2 this line uses element 2 of the y and y1 variables.
k2=f(t+h/2,y(i)+k1*h/2,y1(i)+k1*h/2);
k3=f(t+h/2,y(i)+k2*h/2,y1(i)+k2*h/2);
k4=f(t+h,y(i)+k3*h,y1(i)+k3*h);
y(i+1)=y(i)+(1/6)*(k1+2*k2+2*k3+k4)*h;
At iteration 1 this line assigns to element 2 of y, so at the start of iteration 2 y will have 2 elements.
Nowhere in this loop do you assign to any element of y1, so it will always remain 1 element long. This is the cause of the error at iteration 2. It doesn't make sense to ask for the value of the second element of a vector with one element.
end
If your goal is to solve a system of ODEs, use one of the ones included in MATLAB rather than writing your own.
If your goal is to implement your own ODE solver (because it's a homework assignment, for example) you need to determine what values to assign to subsequent elements in y1 inside the loop like you're assigning to y.
  5 个评论
Steven Lord
Steven Lord 2022-12-15
the y(2) and y(2) value will be gained from the end of the loop of first iteration.
You assign to y(2) in your loop [technically y(i+1) where i is equal to 1.] You never assign to y1(2).
Syed Abdul Rafay
Syed Abdul Rafay 2022-12-15
编辑:Syed Abdul Rafay 2022-12-15
I corrected my code thanks for the help.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Loops and Conditional Statements 的更多信息

产品


版本

R2022b

Community Treasure Hunt

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

Start Hunting!

Translated by