RK4 Errors in MATLAB

2 次查看(过去 30 天)
Kyle Broder
Kyle Broder 2016-8-4
评论: Torsten 2016-8-4

I'm trying to solve the following ODE using RK4 $$F = \frac{1}{50t^2}$$.

My code is the following

    function [t,y]=euler2(y0,dt,t0,tf)
    t=t0:dt:tf;
    y=zeros(1,length(t));
    y(:,1)=y0;
    F = @(t,r) 1/(50*Power(t,2));
    for i = 1:length(t)-1
    k1 = F(y(:,i),t(i));
    k2 = F(y(:,i)+0.5*dt*k1,t(i)+0.5*dt);
    k3 = F(y(:,i)+0.5*dt*k2,t(i)+0.5*dt);
    k4 = F(y(:,i)+dt*k3,t(i)+dt);
    y(:,i+1) = y(:,i) + dt*(k1+2.0*(k2+k3)+k4)/6.0;
    end
    end

When I run the code however, with y0 = 50 for example, or any other number, this yields just an array of zeros and eventually NaN.

Please help.

回答(1 个)

Walter Roberson
Walter Roberson 2016-8-4
You use
F = @(t,r) 1/(50*Power(t,2));
Power appears to be undefined in the code you gave. Perhaps you are thinking of power
Your F function takes t (presumably time) as the first variable and r (unknown purpose) as its second column, and it ignores the second parameter. But you are invoking
k1 = F(y(:,i),t(i));
in which you appear to be passing y as the first parameter and a time as the second parameter. That k1 y(:,i) is going to be received into the parameter that F knows as t, and that k1 t(i) is going to be received into the parameter that F knows as r (where it is going to be ignored).
For the sake of the readers, you need to be consistent about whether it is your first parameter or the second parameter of F that is the time parameter.
  1 个评论
Torsten
Torsten 2016-8-4
... and t0 should be greater than 0.
Best wishes
Torsten.

请先登录,再进行评论。

类别

Help CenterFile Exchange 中查找有关 Programming 的更多信息

Community Treasure Hunt

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

Start Hunting!

Translated by