I am not getting graph for RK4

clear t %Clear old steps

clear y %Clear y values from previous runs

t1=0; % Boundary value 1

t2=10;% Boundary value 2

h=0.1

N=(t2-t1)/h; % Number of steps

y10=1;

y20=0;

y30=0;

%y0=100; % Boundary value of y(a)

%Incremental step size

t(1) = t1; % assign initial location

Y1(1) = y10;

Y2(1) = y20;

Y3(1) = y30;% assign boundary condition

for n=1:N % For loop, sets next t,y values

t(n+1) = t(n)+h;

k1 = fgas(t(n),y1(n)); %Calls the function f(t,y) = dy/dt

k2 = fgas(t(n)+(h/2),y1(n)+(h/2)*k1);

k3 = fgas(t(n)+(h/2),y1(n)+(h/2)*k2);

k4 = fgas(t(n)+(h/2),y1(n)+(h/2)*k3);

y1(n+1)=y1(n)+h*((k1/6)+(k2/3)+(k3/3)+(k4/6))

end

plot(t,Y1)

hold on

%title(['RK4 N=',num2str(N)])

[KSSV]

In the first place is code running?

[Parth Shah]

yes

[KSSV]

No it will show error.....you have to repalce Y with y. To help more, you need to show us what is function fgas.

[Parth Shah]

function fgas=fgas(t,y10)

fgas=-1*y10*y10;

[Parth Shah]

yes

[Parth Shah]

you were correct

[Parth Shah]

can you tell me one more thing?

[KSSV]

What?

[Parth Shah]

clear t %Clear old steps

clear y %Clear y values from previous runs

t1=0; % Boundary value 1

t2=10;% Boundary value 2

h=0.1

N=(t2-t1)/h; % Number of steps

y10=1;

y20=0;

y30=0;

%y0=100; % Boundary value of y(a)

%Incremental step size

t(1) = t1; % assign initial location

Y1(1) = y10;

Y2(1) = y20;

%Y3(1) = y30;% assign boundary condition

for n=1:N % For loop, sets next t,y values

t(n+1) = t(n)+h;

k1 = fgas(t(n),y10(n)); %Calls the function f(t,y) = dy/dt

k2 = fgas(t(n)+(h/2),y10(n)+(h/2)*k1);

k3 = fgas(t(n)+(h/2),y10(n)+(h/2)*k2);

k4 = fgas(t(n)+(h/2),y10(n)+(h/2)*k3);

y10(n+1)=y10(n)+h*((k1/6)+(k2/3)+(k3/3)+(k4/6))

n=1:N % For loop, sets next t,y values

t(n+1) = t(n)+h;

k1 = fcra(t(n),y20(n)); %Calls the function f(t,y) = dy/dt

k2 = fcra(t(n)+(h/2),y20(n)+(h/2)*k10);

k3 = fcra(t(n)+(h/2),y20(n)+(h/2)*k20);

k4 = fcra(t(n)+(h/2),y20(n)+(h/2)*k30);

y20(n+1)=y20(n)+h*((k10/6)+(k20/3)+(k30/3)+(k40/6))

end

plot(t,y10)

hold on

plot(t,Y2)

hold on

i am adding another function

and getting error in the n=1:N % For loop, sets next t,y values

t(n+1) = t(n)+h; line

before fcra

Answer 1

[Alan Stevens]

Looks like you are confusing Y with y.

[Parth Shah]

did not get you

[Alan Stevens]

In these lines

 k1 = fgas(t(n),y1(n)); %Calls the function f(t,y) = dy/dt
    k2 = fgas(t(n)+(h/2),y1(n)+(h/2)*k1);
    k3 = fgas(t(n)+(h/2),y1(n)+(h/2)*k2);
    k4 = fgas(t(n)+(h/2),y1(n)+(h/2)*k3);
    y1(n+1)=y1(n)+h*((k1/6)+(k2/3)+(k3/3)+(k4/6))
   

You have lower case y1. I suspect they should be upper case Y1.

[Parth Shah]

clear t %Clear old steps

clear y %Clear y values from previous runs

t1=0; % Boundary value 1

t2=10;% Boundary value 2

h=0.1

N=(t2-t1)/h; % Number of steps

y10=1;

y20=0;

y30=0;

%y0=100; % Boundary value of y(a)

%Incremental step size

t(1) = t1; % assign initial location

Y1(1) = y10;

Y2(1) = y20;

%Y3(1) = y30;% assign boundary condition

for n=1:N % For loop, sets next t,y values

t(n+1) = t(n)+h;

k1 = fgas(t(n),y10(n)); %Calls the function f(t,y) = dy/dt

k2 = fgas(t(n)+(h/2),y10(n)+(h/2)*k1);

k3 = fgas(t(n)+(h/2),y10(n)+(h/2)*k2);

k4 = fgas(t(n)+(h/2),y10(n)+(h/2)*k3);

y10(n+1)=y10(n)+h*((k1/6)+(k2/3)+(k3/3)+(k4/6))

n=1:N % For loop, sets next t,y values

t(n+1) = t(n)+h;

k1 = fcra(t(n),y20(n)); %Calls the function f(t,y) = dy/dt

k2 = fcra(t(n)+(h/2),y20(n)+(h/2)*k10);

k3 = fcra(t(n)+(h/2),y20(n)+(h/2)*k20);

k4 = fcra(t(n)+(h/2),y20(n)+(h/2)*k30);

y20(n+1)=y20(n)+h*((k10/6)+(k20/3)+(k30/3)+(k40/6))

end

plot(t,y10)

hold on

plot(t,Y2)

hold on

i want to add another function

and i am getting error for second function

[Parth Shah]

why is that so?

[Alan Stevens]

You need to show what function fcra is.

[Parth Shah]

function fcra=fcra(t1,y10,y20)

fcra=((-1*y10*y10)-(0.5*y20));

[Parth Shah]

Index exceeds the number of array elements (2).

Error in RungeKutta2 (line 26)

t(n+1) = t(n)+h;

Answer 2

[Parth Shah]

function fcra=fcra(t1,y10,y20)

fcra=((-1*y10*y10)-(0.5*y20));

[Parth Shah]

Index exceeds the number of array elements (2).

Error in RungeKutta2 (line 26)

t(n+1) = t(n)+h;

[Alan Stevens]

You have given functio fcra three arguments, but you are only calling it with two.

And you don't need to redefine n = 1:N and t(n+1) within your first for n =1:N loop.

[Parth Shah]

okay so what should i write for k1?

[Parth Shah]

k1 = fcra(t(n),y20(n),y10(n)); %Calls the function f(t,y) = dy/dt

k2 = fcra(t(n)+(h/2),y20(n)+(h/2)*k10),y10(n)+(h/2)*k1)

[Parth Shah]

like this way?

[Alan Stevens]

The following works, but you'll need to check carefully that it does what you want:

t1=0;     % Boundary value 1
t2=10;% Boundary value 2
h=0.1;
N=(t2-t1)/h;    % Number of steps
y10=1;
y20=0;
y30=0;
%y0=100;    % Boundary value of y(a)
 %Incremental step size
t(1) = t1;  % assign initial location
Y1(1) = y10;
Y2(1) = y20;
%Y3(1) = y30;% assign boundary condition
for n=1:N  % For loop, sets next t,y values
    t(n+1) = t(n)+h;
   k1 = fgas(t(n),y10(n)); %Calls the function f(t,y) = dy/dt
    k2 = fgas(t(n)+(h/2),y10(n)+(h/2)*k1);
    k3 = fgas(t(n)+(h/2),y10(n)+(h/2)*k2);
    k4 = fgas(t(n)+(h/2),y10(n)+(h/2)*k3);
    y10(n+1)=y10(n)+h*((k1/6)+(k2/3)+(k3/3)+(k4/6));
   
   %n=1:N  % For loop, sets next t,y values
    %t(n+1) = t(n)+h;
    y10a = y10(n);
    y10b = y10(n)+(h/2)*k1; 
    y10c = y10(n)+(h/2)*k2;
    y10d = y10(n)+(h/2)*k3;
    k10 = fcra(t(n),y20(n),y10a); %Calls the function f(t,y) = dy/dt
    k20 = fcra(t(n)+(h/2),y20(n)+(h/2)*k10,y10b);
    k30 = fcra(t(n)+(h/2),y20(n)+(h/2)*k20,y10c);
    k40 = fcra(t(n)+(h/2),y20(n)+(h/2)*k30,y10d);
   y20(n+1)=y20(n)+h*((k10/6)+(k20/3)+(k30/3)+(k40/6));
%      
end
plot(t,y10)
hold on
plot(t,y20)

function fg = fgas(~,y10)
         fg = -1*y10*y10;
end
function fc=fcra(~,y20,y10)
        fc=((-1*y10*y10)-(0.5*y20));
end

[Parth Shah]

graph is very weird

[Alan Stevens]

Possibly so, but it looks like you are trying to solve the following equations:

% dy1/dt = -y1^2    y1(0) = 1
% dy2/dt = -y1^2 - 0.5*y2   y2(0) = 0

If so, using Matlab's ode45 solver we can do the following:

y0 = [1 0];
tspan = [0 10];

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

function dydt = odefn(~,y)

         dydt = [-y(1)^2;
                -y(1)^2 - 0.5*y(2)];
         
end

This results in the same graph:

If you are trying to solve a different set of equations then your code isn't yet finished!

[Parth Shah]

I WNAT THIS GRAPH BUT I WANTED TO DO BY RK4 METHOD

[Alan Stevens]

That's exactly what your code does then!