Runga Kutta method 4

조회 수: 2 (최근 30일)
Ali
Ali 2020년 8월 20일
댓글: Ali 2020년 11월 13일
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
Ali
Ali 2020년 8월 20일
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
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);

댓글을 달려면 로그인하십시오.

채택된 답변

Alan Stevens
Alan Stevens 2020년 8월 20일
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)
  댓글 수: 6
Alan Stevens
Alan Stevens 2020년 11월 13일
Try it. It works, it's just a boring result!
Ali
Ali 2020년 11월 13일
Thanks Alan

댓글을 달려면 로그인하십시오.

추가 답변 (0개)

카테고리

Help CenterFile Exchange에서 Numerical Integration and Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by