CODEWITHSUNDEEP
algorithmmatlabmathnumerical-methodsdifferential-equations
mrc
mrc

Reputation: 3173

How to use the RK4 algorithm to solve an ODE?

I am using an RK4 algorithm:

function R=RK4_h(f,a,b,ya,h)

%  Input
%   - f field of the edo y'=f(t,y). A string of characters 'f'
%   - a and b initial and final time
%   - ya initial value y0
%   - h lenght of the step

%  Output
%   - R=[T' Y'] where T independent variable and Y dependent variable


N = fix((b-a) / h);
T = zeros(1,N+1);
Y = zeros(1,N+1);

%   Vector of the time values 
T = a:h:b;

%   Solving ordinary differential equation
Y(1) = ya;
for j = 1:N
    k1 = h*feval(f,T(j),Y(j));
    k2 = h*feval(f,T(j)+h/2,Y(j)+k1/2);
    k3 = h*feval(f,T(j)+h/2,Y(j)+k2/2);
    k4 = h*feval(f,T(j)+h,Y(j)+k3);
    Y(j+1) = Y(j) + (k1+2*k2+2*k3+k4)/6;
end
R=[T' Y'];

In my main script I call it for every value as:

xlabel('x')
ylabel('y')

h=0.05;
fprintf ('\n First block \n');
xx = [0:h:1];
Nodes = length(xx);
yy = zeros(1,Nodes);
for i=1:Nodes
  fp(i)=feval('edo',-1,xx(i));
end
E=RK4_h('edo',0,1,-1,h);
plot(E);
fprintf ('\n%f',E);

The problem is when I try to use RK4 algorithm with edo formula:

function edo = edo(y,t)
edo = 6*((exp(1))^(6*t))*(y-(2*t))^2+2;

The results are not logical, for example, the real value for are: y(0)=8, y(1)=11,53. But the estimate is not close. Any of both coordinates in E vector represents a feasible approach for the problem, so I do not know if this is the correct implementation.

There is a basic error of implementation?

Upvotes: 0

Views: 84

Answers (1)

Aziz
Aziz

Reputation: 20745

The function edo takes t as the first parameter, and y as the second parameter. You have the parameters reversed.

Your function should be:

 function edo = edo(t,y)   % NOT edo(y,t)
   edo = 6*((exp(1))^(6*t))*(y-(2*t))^2+2;

Upvotes: 1

Related Questions