RK4 Function Code

Question

John Brown 2022-12-5

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1871657-rk4-function-code

编辑： Walter Roberson 2023-8-5

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 个评论
显示无隐藏无

Torsten 2022-12-5

Most probably the reason is that you didn't call "rk4" correctly.

John Brown 2022-12-5

It shouldn’t have a th when called. When testing it didn’t. I forgot to delete that before posting here

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

David Hill 2022-12-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1871657-rk4-function-code#answer_1120627

在 MATLAB Online 中打开

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

2 个评论
显示无隐藏无

John Brown 2022-12-5

For this code i get the error "Error using indexing" and " error in line 6

David Hill 2022-12-5

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?

请先登录，再进行评论。

Answer 2

Steve Areola 2023-8-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1871657-rk4-function-code#answer_1283842

在 MATLAB Online 中打开

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