Skip to main content

FOURTH ORDER RUNGE-KUTTA METHOD- MOTION OF A SPHERICAL MASS WITH AIR RESISTANCE


## Motion of a spherical mass with air resistance
## Fourth order Runge-Kutta Method
## Very Important!!! The positive velocity direction
## is the direction of the gravitational acceleration
m=1E-2; ## mass of the object(kg/m)
g=9.8; ## Acceleration due to the gravity (m/sec^2)
v0=0; ## initial velocity of the object(m/sec)
k=1E-4; ## Air drag coefficient(kg/m)
t=[]; ## Empty time vector(sec)
t(1)=0; ## Released time(sec)
dt=1E-1; ## Increment in time(sec)
v=[]; ## Empty velocity vector(m/s)
v(1)=v0; ## Initial velocity
v_nodrag=[]; ## Velocity by ignoring air drag
## Analytic solution by zeroth order approximation.
v_nodrag(1)=0; ## Velocity vector with no drag.
n=1; ## Initialize the loop index
f0=[];
f1=[]; ##The 1st order derivatives by
f2=[]; ##by Runge-Kutta. Eqn 5.33 pg217
f3=[];
## Run the loop until the time reaches the value 10sec.
while (t(n)<=10);
f0(n)=g-(k/m)*v(n)*v(n);
v_f0(n)=v(n)+(dt/2)*f0(n);
f1(n)=g-(k/m)*v_f0(n)*v_f0(n);
v_f1(n)=v(n)+(dt/2)*f1(n);
f2(n)=g-(k/m)*v_f1(n)*v_f1(n);
v_f2(n)=v(n)+dt*f2(n);
f3(n)=g-(k/m)*v_f2(n)*v_f2(n);
v(n+1)=v(n)+(dt/6)*(f0(n)+2*f1(n)+2*f2(n)+f3(n));
t(n+1)=t(n)+dt;
v_nodrag(n+1)=v_nodrag(n)+g*dt; ##analytic solution
n++;
endwhile
plot(t,v,'r-',t,v_nodrag,'b-');
title('Velocity vs Time');
xlabel('time(sec)');
ylabel('velocity(m/sec)');
legend('v(Runge-Kutta)','v(no drag)')
axis([0,13]);
save -text RUNGEKUTTA.dat
print('-dpsc','RUNGEKUTTA.ps ')

Comments

Post a Comment

Popular posts from this blog

FACTORIAL

## Function that calculates the factorial of a number ## Usage : f=factorial(n) function f=factorial(n) ## Initialize the output f=1; ## Check whether the input is correct if ( (n<0) || (rem(n,1)~=0) ) printf("n cannot be a negative number. Exiting...\n"); return endif for num=1:n f*=num; endfor endfunction
Post a Comment
Read more

One Dimensional Harmonic Oscillator-Numerov Method

x=[]; h0=1; M=4; N=M+1; x(1)=0; x(N)=x(1)+h0*M; x=x(1):h0:x(N) A=zeros(N); A(1,1)=-2*(5*(x(N)*h0)^2/12+1); A(N,N)=-2*(5*(x(1)*h0)^2/12+1); A(1,2)=1-(x(M)*h0)^2)/12; A(N,M)=1-(x(2)*h0)^2)/12; B=zeros(N); B(1,1)=B(N,N)=-10*(h0^2)/6; B(1,2)=B(N,M)=-(h0^2)/6; for i=2:M B(i,i)=-10*(h0^2)/6; B(i,i-1)=B(i,i+1)=-(h0^2)/6; A(i,i)=-2*(5*(x(N+1-i)*h0)^2/12+1); A(i,i+1)=1-(x(N-i)*h0)^2)/12; A(i,i-1)=1-(x(N+2-i)*h0)^2)/12; end A B
Post a Comment
Read more

NEWTON’S METHOD FOR MINIMUM

##Newton's Method to find ##the minimum of the function F(x)=(x-2)^4-9 ##with the initial guess xmin=1.0 ##Constants and initializations xmin=[]; ##The empty array of x that minimizes the F(x) xmin(1)=1.0; ##Initial value of the xmin Fmin=[]; ##Minimum values of F(x) x=0.0:0.1:4.0; ##Only for plotting purposes F=[]; ##Our examined Function evaluated on x-space Fp=[]; ##First derivative of F(x) wrt x Fpp=[]; ##Second derivative o F(x) wrt x NSteps=50; ##Step number of iteration ##Algorithm for n=1:NSteps Fmin(n)=(xmin(n)-2)^4-9; Fp(n)=4*(xmin(n)-2)^3; Fpp(n)=12*(xmin(n)-2)^2; xmin(n+1)=xmin(n)-Fp(n)/Fpp(n); Fmin(n+1)=(xmin(n+1)-2)^4-9; endfor printf("x*, at which F(x) is minimum, is %1.6f\n",xmin(n+1)) printf("Minimum of F(x) is %1.6f\n",Fmin(n+1)) F=(x-2).^4-9; subplot(2,1,1) plot(x,F) title('Newton^,s Method-F(x) vs x'); xlabel('x'); ylabel('F(x)'); text(2,-7,'\downarrow') text(1.7,-5.6,'(xmin,Fm...
Post a Comment
Read more