Runga Kutta method 4
1 vue (au cours des 30 derniers jours)
Afficher commentaires plus anciens
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 commentaires
Réponse acceptée
Alan Stevens
le 20 Août 2020
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 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Numerical Integration and Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!