RK4 Function Code

16 visualizaciones (últimos 30 días)
John Brown
John Brown el 5 de Dic. de 2022
Editada: Walter Roberson el 5 de Ag. de 2023
The issue occurs when running the function at the 'k_1 =' part of the code, the error message says "symbolic function expected 1 input but recieved 2"
function [x,y] = rk4(dydx,xo,xf,yo,h)
x = xo:h:xf ;
y = zeros(1,length(x));
y(1)= yo ;
for i = 1:(length(x)-1)
k_1 = dydx(x(i),y(i));
k_2 = dydx(x(i)+0.5*h,y(i)+0.5*h*k_1);
k_3 = dydx((x(i)+0.5*h),(y(i)+0.5*h*k_2));
k_4 = dydx((x(i)+h),(y(i)+k_3*h));
y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end
dydx = @(x,y) 3.*exp(-x)-0.4*y;
%[x,y] = rk4(dydx,0,100,-0.5,0.5); %plot(x,y,'o-');
end
  2 comentarios
Torsten
Torsten el 5 de Dic. de 2022
Most probably the reason is that you didn't call "rk4" correctly.
John Brown
John Brown el 5 de Dic. de 2022
It shouldn’t have a th when called. When testing it didn’t. I forgot to delete that before posting here

Iniciar sesión para comentar.

Respuestas (2)

David Hill
David Hill el 5 de Dic. de 2022
dydx = @(x,y) 3.*exp(-x)-0.4*y;
[x,y] = rk4(dydx,0,100,-0.5,0.5);
plot(x,y,'o-');
function [x,y] = rk4(dydx,xo,xf,yo,h)
x = xo:h:xf ;
y = zeros(1,length(x));
y(1)= yo ;
for i = 1:(length(x)-1)
k_1 = dydx(x(i),y(i));
k_2 = dydx(x(i)+0.5*h,y(i)+0.5*h*k_1);
k_3 = dydx((x(i)+0.5*h),(y(i)+0.5*h*k_2));
k_4 = dydx((x(i)+h),(y(i)+k_3*h));
y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;
end
end
  2 comentarios
John Brown
John Brown el 5 de Dic. de 2022
For this code i get the error "Error using indexing" and " error in line 6
David Hill
David Hill el 5 de Dic. de 2022
The above code runs without failure as you can see above. Did you save the function and run the 3-line script in the command line?

Iniciar sesión para comentar.


Steve Areola
Steve Areola el 5 de Ag. de 2023
clc
a=0;b=2;N=10;alph=0.5;
h=(b-a)/N;t(1)=a;w(1)=alph;
f=@(t,y) y-t^2+1;
for i=2:N+1
k1= h*f(t(i-1),w(i-1));
k2= h*f(t(i-1)+(h/2),w(i-1)+(k1/2));
k3= h*f(t(i-1)+(h/2),w(i-1)+(k2/2));
k4= h*f(t(i-1)+h,w(i-1)+k3);
w(i)=w(i-1)+(k1+2*k2+2*k3+k4)/6;
t(i)=a+(i-1)*h;
end
x=[t;w];
disp("t w")
t w
fprintf("%5.7f %5.7f\n",x)
0.0000000 0.5000000 0.2000000 0.8292933 0.4000000 1.2140762 0.6000000 1.6489220 0.8000000 2.1272027 1.0000000 2.6408227 1.2000000 3.1798942 1.4000000 3.7323401 1.6000000 4.2834095 1.8000000 4.8150857 2.0000000 5.3053630

Categorías

Más información sobre Programming en Help Center y File Exchange.

Etiquetas

Community Treasure Hunt

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

Start Hunting!

Translated by