Runga Kutta method 4

Question

Ali 2020-8-20

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/581963-runga-kutta-method-4

评论： Ali 2020-11-13

采纳的回答： Alan Stevens

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

3 个评论
显示 1更早的评论隐藏 1更早的评论

Ali 2020-8-20

coupled Equation by RK method.m

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 2020-11-13

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 2020-8-20

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/581963-runga-kutta-method-4#answer_482549

在 MATLAB Online 中打开

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)