Solving differential equation with Runge Kutta 4th order

I've tried to write this code but it seems to give different values than ODE 45.

The equation is the following:

Does it make sense?

Thanks in advance

%Extremes

a=0;

b=1;

%Stepsize

h=0.05;

%time interval and initial value

t=(a:h:b)';

x(1)=0.032;

parameter=0.15;

t1=zeros(size(t));

n=numel(t1);

%RK cycle

for i=1:n-1

%function

dxdt = @(x) parameter*(1-x)-(1-x)*(35*x(i)/((x(i)^2)/0.01+x(i)+1));

j1 = h*dxdt(x(i));

j2 = h*dxdt(x(i)+j1/2);

j3 = h*dxdt(x(i)+j2/2);

j4 = h*dxdt(x(i)+j3);

x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4);

end

Answer 1

[VBBV]

a=0;
b=1;
%Stepsize
h=0.05;
%time interval and initial value
t=(a:h:b)';
x(1)=0.032;
parameter=0.15;
t1=zeros(size(t)); 
n=numel(t1);
%RK cycle
dxdt = @(x) parameter*(1-x)-(1-x).*(35*x./((x.^2)/0.01+x+1));
for i=1:n-1 
%function
j1 = h*dxdt(x(i));
j2 = h*dxdt(x(i)+j1/2); 
j3 = h*dxdt(x(i)+j2/2);
j4 = h*dxdt(x(i)+j3);
x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4); 
end
plot(t,x)

[Edoardo Bertolotti]

Thank you. I should have test my code before posting it, I'm modifying a more complex problem with more than 1 equation (but they are all indipendent).

Now I was comparing the methods (also with Euler) and there's still difference. I've tried to simplify the problem again but now ODE45 does not work. Yet, in the original problem ODE45 gives the same results as Euler..so I still think there's something which is not correct.

Here I report the simplified version where ODE45 does not work but Euler yes.

a=0;

b=6;

%Stepsize

h=0.05;

%time interval and initial value

t=(a:h:b)';

x(1)=0.3075;

parameter=1.5;

t1=zeros(size(t));

n=numel(t1);

%Euler

dxdt = @(x) parameter*(1-x)-(1-x).*(35*x./((x.^2)/0.01+x+1));

for i=1:n-1

x(i+1)=x(i)+h*dxdt(x(i));

end

ylim([-0.1 1.1])

hold on

plot(t,x)

title('methods')

hold on

grid on

xlabel('time')

ylabel('concentration')

hold on

%RK cycle

dxdt = @(x) parameter*(1-x)-(1-x).*(35*x./((x.^2)/0.01+x+1));

for i=1:n-1

%function

j1 = h*dxdt(x(i));

j2 = h*dxdt(x(i)+j1/2);

j3 = h*dxdt(x(i)+j2/2);

j4 = h*dxdt(x(i)+j3);

x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4);

end

hold on

plot(t,x)

hold on

run trial_ODE45_solver.m

legend('Euler','Runge Kutta','ODE45','Location','NorthOutside','Orientation','horizontal','Box','off')

Here's trial_ODE45_solver.m

y0 = 0.3075;

tspan = [0 6];

[t,y] = ode45(@trial_ODE45,tspan,y0);

plot(t,y,'-')

and trial_ODE45

function dydt = ODEsystem_3(t,y)

parameter = 1.5;

%equations system

dydt = @(y) parameter*(1-y)-(1-y).*(35*y./((y.^2)/0.01+y+1));

end

[Rakesh]

Hi

In the above solution using rk4 method. I want to give jump to the solution at discrete points then how will I code this?

Many thanks in advance!

[Torsten]

Integrate to the first jump point,

set the new initial values to the solution at the last time instant + the jump in the solution,

restart and integrate to the next jump point,

set the new initial values to the solution at the last time instant + the jump in the solution,

...

[Ahmed J. Abougarair]

clc;

clear all;

F = @(t,y) 4*exp(0.8*t)-0.5*y

t0=input('Enter the value of t0 : \n');

y0=input('Enter the value of y0 : \n');

tn=input('Enter the value of t for which you want to find the value of y : \n');

h=input('Enter the step length : \n');

i=0;

while i<tn

k_1 = F(t0,y0);

k_2 = F(t0+0.5*h,y0+0.5*h*k_1);

k_3 = F((t0+0.5*h),(y0+0.5*h*k_2));

k_4 = F(((t0)+h),(y0+k_3*h));

nexty = y0 + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h;

y0=nexty;

t0=t0+h;

i=i+h;

end

fprintf('The value of y at t=%f is %f \n',t0,y0)

% validate using a decent ODE integrator

tspan = [0,1]; Y0 = 2;

[tx,yx] = ode45(F, tspan, Y0);

fprintf('The true value of y at t=%f is %f \n',tspan(end),yx(end))

Et= (abs(yx(end)-y0)/yx(end))*100;

fprintf('The value of error Et at t=%f is %f%% \n',tspan(end),Et)

Answer 2

[Torsten]

y0 = 0.3075; 
tspan = [0 6]; 
[t,y] = ode45(@trial_ODE45,tspan,y0);
plot(t,y,'-')

function dydt = trial_ODE45(t,y)
parameter = 1.5;
%equations system  
dydt = parameter*(1-y)-(1-y).*(35*y./((y.^2)/0.01+y+1));
end

Answer 3

[Edoardo Bertolotti]

Ok, thanks to both, the code is working now, but still, like the main problem, it shows the same results:

Euler is similar to ODE45, but RK4 no..I would expect ODE45 to be more similar to RK than Euler or equal to both, using same steps etc

a=0;

b=6;

%Stepsize

h=0.05;

%time interval and initial value

t=(a:h:b)';

x(1)=0.3075;

parameter=1.5;

t1=zeros(size(t));

n=numel(t1);

%Euler

dxdt = @(x) parameter*(1-x)-(1-x).*(35*x./((x.^2)/0.01+x+1));

for i=1:n-1

x(i+1)=x(i)+h*dxdt(x(i));

end

ylim([-0.1 1.1])

hold on

plot(t,x)

title('methods')

hold on

grid on

xlabel('time')

ylabel('concentration')

hold on

%RK cycle

dxdt = @(x) parameter*(1-x)-(1-x).*(35*x./((x.^2)/0.01+x+1));

for i=1:n-1

%function

j1 = h*dxdt(x(i));

j2 = h*dxdt(x(i)+j1/2);

j3 = h*dxdt(x(i)+j2/2);

j4 = h*dxdt(x(i)+j3);

x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4);

end

hold on

plot(t,x)

hold on

%ODE45

y0 = 0.3075;

tspan = [0 6];

[t,y] = ode45(@trial_ODE45,tspan,y0);

plot(t,y,'-')

legend('Euler','Runge Kutta','ODE45','Location','NorthOutside','Orientation','horizontal','Box','off')

% % %runge kutta

function dydt = trial_ODE45(t,y)

parameter = 1.5;

%equations system

dydt = parameter*(1-y)-(1-y).*(35*y./((y.^2)/0.01+y+1));

end

[Torsten]

x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3+j4);

instead of

x(i+1) = x(i)+(1/6)*(j1+2*j2+2*j3*j4);

in RK4.

[Edoardo Bertolotti]

Oh, amazing, thank you!