Runga Kutta method 4

I want to solve this equation by Runga kutta method 4.

d2xd/dt2+0.2dx/dt-2.45xd+0.2908xp=

with intial condition

t=to=0

xd=1

xp=0

from t=0 to 100

[James Tursa]

What have you done so far? What specific problems are you having with your code? Are you allowed to use ode45, or are you required to write your own RK4 code?

[Ali]

This is my code but it does not run.This code is for two equation of similar kind.But if one equation code run that s also ok.

[Ali]

I want to solve second order equation through Ruga kutta method .Please check what is error in this code

m=100;

k=25000;

e=100;

F=15000;

w=(k/m)^0.5;

t(1)=0;

x(1)=0;

v(1)=1;

h=0.1;

tfinal=100;

N=ceil((tfinal-t(1))/h);

f1=@(t,x,v) v;

f2=@(t,x,v) -(k/m)*v-(e/m)*x+(F/m)*cos(w*t);

for j=1:N;

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

K1x=f1(t(j),x(j),v(j));

K1v=f2(t(j),x(j),v(j));

K2x=f1(t(j)+h/2,x(j)+h/2*K1x,v(j)+h/2*K1v);

K2v=f2(t(j)+h/2,x(j)+h/2*K1x,v(j)+h/2*K1v);

K3x=f1(t(j)+h/2,x(j)+h/2*K2x,v(j)+h/2*K2v);

K3v=f2(t(j)+h/2,x(j)+h/2*K2x,v(j)+h/2*K2v);

K4x=f1(t(j)+h,x(j)+h*K3x,v(j)+h*K3v);

K4v=f2(t(j)+h,x(j)+h*K3x,v(j)+h*K3v);

x(j+1)=x(j)+h/6*(K1x+2*K2x+2*K3x+K4x);

v(j+1)=v(j)+h/6*(K1v+2*K2v+2*K3v+K4v);

end

disp(x(N));

figure;

plot(t,x,'lineWidth',1);

xlabel('t');

ylabel('X','V');

hold on

vel(v(N));

figure;

plot(t,v,lineWidth',1);

Answer 1

[Alan Stevens]

Firstly, you can't have spaces in the name of the m-file. Remove them.

Secondly, remember "length" not "lenght".

Thirdly, pass parameters to the functions in the same order in which you define the functions.

Fourthly, make use of Matlab's error messages!

Here is a working code, but it produces nonsense because of the third point above. Put the parameters in the right order and you might get something sensible.

%input paramters
h=1;
t=1:h:100;
xd=zeros(1,length(t));  %%%%% length not lenght
xp=zeros(1,length(t));  %%%%% length not lenght
vd=zeros(1,length(t));  %%%%% length not lenght
vp=zeros(1,length(t));  %%%%% length not lenght
xdo=1;
vdo=0;
xpo=-1;
vpo=0;
xd(1)=xdo; %%%%% These next four lines must come after the previous four, 
vd(1)=vdo; %%%%% and array indices start at 1 not 0 in Matlab.
xp(1)=xpo;
vp(1)=vpo;


%RK4
f1=@(t,xd,xp,vd,vp) vd;
f2=@(t,xd,xp,vd,vp) -0.4198*vd-5.9282*xd-1.5285*xp;
f3=@(t,xd,xp,vd,vp) vp;
f4=@(t,xd,xp,vd,vp) -0.4047*vp+1.5285*xp-7.3774*xd; %%%% xd not xp ?


for i=1:(length(t)-1)   %%%%% length not lenght
    
   %%%% The following parameters are not passed to the functions f1, f2,
   %%%% etc in the order in which they are defined.  You need to correct
   %%%% that.
    K1xd=f1(t(i),xd(i),vd(i), xp(i),vp(i));
    K1vd=f2(t(i),xd(i),vd(i),xp(i),vp(i));
    K1xp=f1(t(i),xp(i),vp(i), xd(i),vd(i));
    K1vp=f2(t(i),xp(i),vp(i),xd(i),vd(i));
    
    K2xd=f1(t(i)+h/2,xd(i)+h*K1xd/2,vd(i)+h*K1vd/2,xp(i)+h*K1xp/2,vp(i)+h*K1vp/2);
    K2vd=f2(t(i)+h/2,xd(i)+h*K1xd/2,vd(i)+h*K1vd/2,xp(i)+h*K1xp/2,vp(i)+h*K1vp/2); %%%%%%%%%
    K2xp=f1(t(i)+h/2,xp(i)+h*K1xp/2,vp(i)+h*K1vp/2,xd(i)+h*K1xd/2,vd(i)+h*K1vd/2);
    K2vp=f2(t(i)+h/2,xp(i)+h*K1xp/2,vp(i)+h*K1vp/2,xd(i)+h*K1xp/2,vd(i)+h*K1vd/2); %%%%%%%%%
 
    K3xd=f1(t(i)+h/2,xd(i)+h*K2xd/2,vd(i)+h*K2vd/2,xp(i)+h*K2xp/2,vp(i)+h*K2vp/2);
    K3vd=f2(t(i)+h/2,xd(i)+h*K2xd/2,vd(i)+h*K2vd/2,xp(i)+h*K2xp/2,vp(i)+h*K2vp/2); 
    K3xp=f1(t(i)+h/2,xp(i)+h*K2xp/2,vp(i)+h*K2vp/2,xd(i)+h*K2xd/2,vd(i)+h*K2vd/2);
    K3vp=f2(t(i)+h/2,xp(i)+h*K2xp/2,vp(i)+h*K2vp/2,xd(i)+h*K2xd/2,vd(i)+h*K2vd/2); 
    
    K4xd=f1(t(i)+h,xd(i)+h*K3xd,vd(i)+h*K3vd,xp(i)+h*K3xp,vp(i)+h*K3vp);
    K4vd=f2(t(i)+h,xd(i)+h*K3xd,vd(i)+h*K3vd,xp(i)+h*K3xp,vp(i)+h*K3vp);   
    K4xp=f1(t(i)+h,xp(i)+h*K3xp,vp(i)+h*K3vp,xd(i)+h*K3xd,vd(i)+h*K3vd);
    K4vp=f2(t(i)+h,xp(i)+h*K3xp,vp(i)+h*K3vp,xd(i)+h*K3xd,vd(i)+h*K3vd);    %%%%%%%  K4vp not K4vd
    
    xd(i+1)=xd(i)+h/6*(K1xd+K1xp+2*K2xd+2*K2xp+2*K3xd+2*K3xp+K4xd+K4xp);
    vd(i+1)=vd(i)+h/6*(K1vd+K1vp+2*K2vd+2*K2vp+2*K3vd+2*K3vp+K4vd+K4vp);
    xp(i+1)=xp(i)+h/6*(K1xp+K1xd+2*K2xp+2*K2xd+2*K3xp+2*K3xd+K4xp+K4xd);
    vp(i+1)=vp(i)+h/6*(K1vp+K1vd+2*K2vp+2*K2vd+2*K3vp+2*K3vd+K4vp+K4vd);
end
%plots
plot(t,xd)
hold on
plot(t,xp)

[Ali]

Thanks

[Ali]

I want to solve second order equation through Ruga kutta method .Please check what is error in this code

m=100;

k=25000;

e=100;

F=15000;

w=(k/m)^0.5;

t(1)=0;

x(1)=0;

v(1)=1;

h=0.1;

tfinal=100;

N=ceil((tfinal-t(1))/h);

f1=@(t,x,v) v;

f2=@(t,x,v) -(k/m)*v-(e/m)*x+(F/m)*cos(w*t);

for j=1:N;

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

K1x=f1(t(j),x(j),v(j));

K1v=f2(t(j),x(j),v(j));

K2x=f1(t(j)+h/2,x(j)+h/2*K1x,v(j)+h/2*K1v);

K2v=f2(t(j)+h/2,x(j)+h/2*K1x,v(j)+h/2*K1v);

K3x=f1(t(j)+h/2,x(j)+h/2*K2x,v(j)+h/2*K2v);

K3v=f2(t(j)+h/2,x(j)+h/2*K2x,v(j)+h/2*K2v);

K4x=f1(t(j)+h,x(j)+h*K3x,v(j)+h*K3v);

K4v=f2(t(j)+h,x(j)+h*K3x,v(j)+h*K3v);

x(j+1)=x(j)+h/6*(K1x+2*K2x+2*K3x+K4x);

v(j+1)=v(j)+h/6*(K1v+2*K2v+2*K3v+K4v);

end

disp(x(N));

figure;

plot(t,x,'lineWidth',1);

xlabel('t');

ylabel('X','V');

hold on

vel(v(N));

figure;

plot(t,v,lineWidth',1);

[Alan Stevens]

Try the following

m=100;
k=25000;
e=100;
F=15000;
w=(k/m)^0.5;
t(1)=0;
x(1)=0;
v(1)=1;
h=0.01;              %%%% Original timetep was too large
tfinal=1;            %%%% Not much point in going as far as 100.
N=ceil((tfinal-t(1))/h);
f1=@(t,x,v) v;
f2=@(t,x,v) -(k/m)*v-(e/m)*x+(F/m)*cos(w*t);
for j=1:N
    t(j+1)=t(j)+h;
    K1x=f1(t(j),x(j),v(j));
    K1v=f2(t(j),x(j),v(j));
    K2x=f1(t(j)+h/2,x(j)+h/2*K1x,v(j)+h/2*K1v);
    K2v=f2(t(j)+h/2,x(j)+h/2*K1x,v(j)+h/2*K1v);
    K3x=f1(t(j)+h/2,x(j)+h/2*K2x,v(j)+h/2*K2v);
    K3v=f2(t(j)+h/2,x(j)+h/2*K2x,v(j)+h/2*K2v);
    K4x=f1(t(j)+h,x(j)+h*K3x,v(j)+h*K3v);
    K4v=f2(t(j)+h,x(j)+h*K3x,v(j)+h*K3v);
    x(j+1)=x(j)+h/6*(K1x+2*K2x+2*K3x+K4x);
    v(j+1)=v(j)+h/6*(K1v+2*K2v+2*K3v+K4v);
end
disp(x(N));
figure;
plot(t,x);
xlabel('t');
ylabel('X, V');
hold on
disp(v(N));
figure;
plot(t,v);

[Ali]

Thank

what changes are made if I increase time upto to 100 sec.

[Alan Stevens]

Try it. It works, it's just a boring result!

[Ali]

Thanks Alan