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

SIMPLE AND MODIFIED EULER METHOD

##Usage:Call Octave from terminal ##and then call EulerMethodUmitAlkus.m ##from octave and finally ##press enter. That's all. ##Simple Euler Method ##Constants and initializations x=[]; ## initial empty vector for x y=[]; ## initial empty vector for y x(1)=1; ## initial value of x y(1)=1; ## initial value of y h=1E-3; ## increment in x dery=[]; ## 1st derivative of y wrt x dery(1)=0;## 1st entry of dery n=1; ## inital loop index for while ## enter the while loop for the interval x=[1,2] while (x(n)<=2) x(n+1)=x(n)+h; dery(n+1)=x(n)*x(n)-2*y(n)/x(n); ##given y(n+1)=y(n)+h*dery(n+1); ##Euler method n++; endwhile ##exit from the 1st while loop ##Modified Euler Method ##Constant and initializations x(1)=1; ## beginnig of the interval [1,2] ymod(1)=1; ## inital value for modified y. ymid=[]; ## empty vector function evaluated at x midpoint xmid=[]; ## empty vector func. of midpoints of the interval h in x-axis. derymod=[]; ## modified derivatives of ymod
Post a Comment
Read more

NEWTON-RAPSON METHOD-8th degree Legendre polynomial

## Newton-Rapson Method to the smallest non negative root ## of the 8th degree Legendre Polynomial ## P8(x)=(1/128)(6435x^8-12012x^6+6930x^4-1260x^2+35) ## where -1<=x<=1. ## for the smallest non negative root, we can ignore ## all the terms except the last two by truncated ## the function to be zero and find ## x=0.167 as the initial smallest non negative ## root. ##Constants and initializations x=[]; ## Empty array for the iterated x roots x(1)=0.16700000; ## Initial guess to begin the iteration for the ## smallest non-negative root. L8=[]; ## Empty array for the Legendre polynomial L8p=[]; ## Empty array for the derivative of the Legendre polynomial for i=1:100 ##The value of the function at x L8(i)=(1/128)*(6435*x(i)^8-12012*x(i)^6+6930*x(i)^4-1260*x(i)^2+35); ##The value of the derivative of the function at x L8p(i)=(1/128)*(6435*8*x(i)^7-12012*6*x(i)^5+6930*4*x(i)^3-1260*2*x(i)); x(i+1)=x(i)-L8(i)/L8p(i); ## the iteration endfor ## For plot let's
Post a Comment
Read more