Skip to main content

Posts

Showing posts from July, 2016

Second Harmonic Generation N=1:21

gnuplot> set xrange [-180:180] gnuplot> set yrange [-180:180] splot sin(cos(x*pi/180))*sin(cos(x* pi/180))/(cos(x*pi/180)*cos(x* pi/180))*sin(cos(y*pi/180))* sin(cos(y*pi/180))/(cos(y*pi/ 180)*cos(y*pi/180))+sin(cos(( x+1)*pi/180))*sin(cos((x+1)* pi/180))/(cos((x+1)*pi/180)* cos((x+1)*pi/180))*sin(cos((y+ 1)*pi/180))*sin(cos((y+1)*pi/ 180))/(cos((y+1)*pi/180)*cos(( y+1)*pi/180))+sin(cos((x+2)* pi/180))*sin(cos((x+2)*pi/180) )/(cos((x+2)*pi/180)*cos((x+2) *pi/180))*sin(cos((y+2)*pi/ 180))*sin(cos((y+2)*pi/180))/( cos((y+2)*pi/180)*cos((y+2)* pi/180))+sin(cos((x+3)*pi/180) )*sin(cos((x+3)*pi/180))/(cos( (x+3)*pi/180)*cos((x+3)*pi/ 180))*sin(cos((y+3)*pi/180))* sin(cos((y+3)*pi/180))/(cos(( y+3)*pi/180)*cos((y+3)*pi/180) )+sin(cos((x+4)*pi/180))*sin( cos((x+4)*pi/180))/(cos((x+4)* pi/180)*cos((x+4)*pi/180))* sin(cos((y+4)*pi/180))*sin( cos((y+4)*pi/180))/(cos((y+4)* pi/180)*cos((y+4)*pi/180)) +sin(cos((x+5)*pi/180))*sin( cos((x+5)*pi/180))/(cos((x+5)* pi/180)*cos((

Second Harmonic Generation

gnuplot> set xrange [-180:180] gnuplot> set yrange [-180:180] gnuplot> set pm3d gnuplot> set hidden3d  gnuplot> set title 'SHG' gnuplot> splot sin(cos(x*pi/180))*sin(cos(x* pi/180))/(cos(x*pi/180)*cos(x* pi/180))*sin(cos(y*pi/180))* sin(cos(y*pi/180))/(cos(y*pi/ 180)*cos(y*pi/180)) title 'N=1 '

FREE FALL

## A function, free fall, that takes in h (in metres), and ## returns the final velocity of the ball at the ## time step just before it touches the ground and makes ## a plot of velocity versus time. function freefall(h) ## constants and initializations v=0; ## initial velocity(at rest) [m/sec] y=h; ## initial altitude of the ball from the ground [m] t=0; ## initial time[sec ] B1=0.05; ## Coeff of the term prop to v [kg/sec] B2=6E-4; ## Coeff of the term prop to v^2 [kg/m] m=0.25; ## Mass of the ball [kg] g=9.8; ## Gravitational acceleration [m/sec^2] dt=0.1; ## Time increment [sec] n=1; ## Initialize the loop index[dimensionless] ## Run the loop or iterate until the vertical component is ## smaller than zero. while (y(n)>0); ## decrease the altitude in each step y=[y;y(n)+v(n)*dt]; ## and accumlate the results in the array ## of the vertical displacement. ## increase the time in each step t=[t;t(n)+dt]; ## and accumlate the results in the array ## of the tim

INTEGRAL (TRAPEZOIDAL RULE)

## This m-file computes the integral of a function f over the ## equally spaced points x using the trapezoidal rule. function I = trapezoid(x,f) h = x(2)-x(1); N = length(f); I = 0; for n=2:N-1, I = I + h*f(n); endfor I = I + h/2*(f(1) + f(N)); endfunction

TOTAL ISING ENERGY

function E=Ising_Energy(lattice,i,j,rneigh) H=2; ## applied magnetic field times permeability[joule] E=0; ## initial configuration energy[joule] J0=1 ## coupling constant of function J, normaly in joule. alpha=2 ## exponential constant of function J. L=length(lattice); ## number of particles or spins in lattice for r=1:rneigh NN1=mod(j+r,L); NN1+=(NN1==0)*L; NN2=mod(j-r,L); NN2+=(NN2==0)*L; NN3=mod(i+r,L); NN3+=(NN3==0)*L; NN4=mod(i-r,L); NN4+=(NN4==0)*L; J=coupling(rneigh,J0,alpha); dE=J*(lattice(i,NN1)+lattice(i,NN2)+ lattice(NN3,j)+lattice(NN4,j))-H*(sum((rand(L)>0.5)*2-1)); E+=-1*dE endfor endfunction

NONUNFORM SOUND WAVE GENERATOR

## a script for a sound wave source made up of two different masses ##seperated in the middle. Let the rigth part is lighter than the ##left one that ensures that the ligther is faster, and so has ##greater r since r=cdt/dx which is dimensionless(ratio). dx=1e-2; ## Increment in the horizontal distance (m) c1=300; ## speed of the left one (m/s) c2=500; ## Speed of the rigth one (m/s) dt=dx/max(c1,c2); ## Increment in time (s). The max c, the min dt. r1=c1*dt/dx; ## Ratio r for the left one r2=c2*dt/dx; ## Ration r for the right one x=-1:dx:1; ## same as the string_fixed l=length(x); x0=0.5; k=1e2; ## Inital profile(return to guassuian distibution) y=initial_profile(x,x0,k); axis([-1.05,1.05,-1.1,1.1]) plot(x,y,'r;;',[0 0],[-1.1 1.1],'r-;;') pause ## Boundary conditions while t=dtn=0 ynow=y; yprev=y; N=100; for n=1:N ynext=propagate_two_parts(ynow,yprev,r1,r2); ## calling the prepared ## function plot(x,ynext,';;',[0 0],[-1.1 1.1],'r

STIFFNESS OF BENDING STRINGS

## A 'realistic', 'non-elastic' string, which responses to any ## bending and has stifness. This script takes in the previous ## and the present profiles and iterates to find ## the profile in the next time step. The ratio 'r' is not 1 ## like in the 'non-realistic' string since the speed of the wave ## always less than the speed of the string it should be less than 1 ## for best and most stable solution ##constants dx=1e-2 ## Spatial increment (m) L=2 ## Length of the string (m) M=L/dx ## Dimensionless partition r=0.25 ## Famous dimensionless ratio E=1e-4 ## Dimensionless stiffnes x=-1:dx:1; l=length(x); x0=0.5; k=1e2; ## Set up the initial profile y=initial_profile(x,x0,k); plot(x,y) pause ## Impose the time boundary condition ynow=y; yprev=y; ynow(1)=ynow(2)=0 Nsteps=2000; for n=3:Nsteps ynow(Nsteps-1)=ynow(Nsteps-2)=0; ynext=propagate_stiff(ynow,yprev,r); plot(x,ynext,';;') axis([-1.05,1.05,-1.1,1.1]) pause(0);

TRAVELING WAVE ON A NONREALISTIC STRING

## Script that simulates a traveling wave on an nonrealistic string. ##The one end of the string is driven sinusoidally ## while the other end is kept fixed. dx=1e-2; ## increment in horizontal displacement(m) c=300; ## speed of the wave (m/s) dt=dx/c; ## Increment in time (s) r=c*dt/dx; ## Ratio of speed of wave to ## the speed of the string ## Initial profile x=-1:dx:1; y=zeros(size(x)); ## Change the gaussian distribution ## to get zero vertical ##displacements for all parts of ##the string.(purely flat initially) plot(x,y) pause ##For the sinusoidally driven end, given constans are: A=0.5; ##amplitude(m) rate=100; ## inverse of the period (cycles/s) omega=2*pi*rate ; ## angular frequency(Radians/s) ## Time boundary conditions. Get previous and ## present displacement as initial profile ynow=y; yprev=y; Nsteps=2000; ## all steps are the same for the loop as the string_fixed ## apart from the function ynext for n=1:Nsteps ynext=propagate_driven(ynow,yprev

BIO-SAVART LAW

/*A C code that returns the magnetic field at a mesh of points *on the xy-plane (xp,yp) for a straight wire segment through which *current passes. The total magnetic field used in this code *has the form (cos(theta1)-cos(theha2))/(yp*r_mag) *ignoring the constants permeability, pi and current I and *taking them together as 1. Thus no need for trapezoid since *the integral of Bio-Savart can be tractable!!! *r_mag represents the distance between the dl segment and *the point (xp,yp) where we will find the B*/ #include <stdio.h> #include <math.h> /*Global constant*/ #define L 3.0 /* length of the wire*/ double r_mag(double x, double xp, double yp) { double rx; rx=xp-x; /* x is any point on xprime*/ return sqrt(rx*rx+yp*yp); } main() { int i; double xp,yp; double B; /*the z component of mag-field that is total field.*/ FILE *fid; fid=fopen("mag-field","w"); /*loop over xp and yp values*/ for (yp=-2.0;yp<=2.0;yp=yp+0.1) for (xp=-

TOTAL DISPLACEMENT OF LOOPLESS RANDOM WALK

## A for good and evil as well as loopless ## total displacement function of random walk ## that takes in the step numbers as an ## argument and returns the displacement ## of the rnd walker.Note that the probabilities ## to turn rigth and left are equal but the left ## step is twice the rigth step. ## Usage : rw_uneven(N) function rw_uneven(N) rn=rand(N,1); ## If the element of a random vector r=-(rn<0.5); ## is smaller than 0.5, then return -1, li=find(r == 0); ## Else return 2(equal probability). r(li)=2; x=sum(r) ## total displacement endfunction

MEAN SQUARED DISPLACEMENT-1D RANDOM WALK

## Routine random walk(rw) for simulating ## one-dimensional random walk ## and calculating the mean squared displacement. x2ave=[]; ## Initial array of the squared ## displacement at initial time(meters). Nrw=3000; ## Number of walkers (Dimensionless). beg=1000; ## Minimum step number ( // ). inc=1000; ## Step number increment ( // ). Nstepsmax=10000; ## Maximum step number ( // ). dt=1; ## Time increment(second). ##Entering the nested loops. for Nsteps=beg:inc:Nstepsmax; ## The loop through the number Nsteps ## of steps each in walk. x2=0; ## Initial squared location of all of the ## walkers at all steps. for m=1:Nrw ## The loop through the desired ## number of walkers. r=rand_disc_rev(Nsteps,0.5); ## Calling the rw generator ## function which will give column vector ## of each elements are either 1 or -1. x2+=sum(r)^2; ## Total displacement squared whose ## elements stems from the vector r. endfor x2ave=[x2ave;x2/Nrw]; ## Accumulation of the squared displacem

RANDOM WALK GENERATOR

## This function ## provides us with very simple ## random walk generator compared ## to the function rand_disc_loop ## which is time-consuming, in operation, ## long, and unuseful for other than the ## coin toss simulatons function xy=rand_disc(N); ## RW generator. r=floor(rand(N,1)*4); ## Random coulumn vector of N ## integer elements ## multiplied by 4 to widen the ## interval from [0,1] to [0,3]. x=y=zeros(size(r)); ## The coulumn vectors ## of N elements with all elements zero. x(find(r==0)) = 1; ## The elements of x, ## the vector function ## which takes in the row # ## of the zero elements of ## the vector r as an argument, ## is assigned to 1. x(find(r==1)) =-1; ## The elements of x, ## the vector function ## which takes in the row # of the ## elements of 1 of the vector r ## as an argument, ## is assigned to -1. y(find(r==2)) = 1; ## The elements of y, ## the vector function ## which takes in the row # of the ## elements of 2 of the vector r ## a

PROPAGATION OF NON-UNIFORM STRING

## The string is made up of two different mass density strings ## seperated in the middle. Choose the right one to be ligther ## than the left. So, since velocity is inversely propotional ## to the mass density, than the rigth one is faster also. ## r1=c1*dt/dx and r2=c1*dt/dx function ynext=propagate_two_parts(ynow,yprev,r1,r2) ## initiallization for, take y as previos one ynext=ynow; ## we have two parts, hence two loops are needed ## the first part is between x=0 and ## let y(now)/2 for better visiulation. ## floor(x) returns the largest integer not greater than x ## length(x) determines the number of column ## or rows in matrix or vector ##I started the 'for' from x=0 but didn't work ##then started from x=1 ## but in this case the ## string_two_parts didn't work ## let x0=i0=2 for i=2:floor(length(ynow)/2) ynext(i) = 2*(1-r1^2)*ynow(i)-yprev(i)+r1^2*(ynow(i+1)+ynow(i-1)); endfor ## the second part is between x0=(L1+L2) /2 and L1+L2 where ## L

REALISTIC AND NON-REALISTIC STRINGS

## A 'realistic', 'non-elastic' string, which responses to any ## bending and has stifness. This script takes in the previous ## and the present profiles and iterates to find ## the profile in the next time step. The ratio 'r' is not 1 ## like in the 'non-realistic' string since the speed of the wave ## always less than the speed of the string it should be less than 1 ## for best and most stable solution ##constants dx=1e-2 ## Spatial increment (m) L=2 ## Length of the string (m) M=L/dx ## Dimensionless partition E=1e-4 ## Dimensionless stiffnes function ynext=propagate_stiff(ynow,yprev,r) ## Quick and dirty way to fix boundary conditions -- for each step ## they are the same as the previous step. ynext=ynow; ynow(1)=ynow(2)=0; ##Entering the loop for i=3:length(ynow)-1 ## boundary condition ynow(length(ynow)-1)=ynow(length(ynow)-2)=0; ## Divide the ynext with many terms into three parts for easiness ynext(i)=(2−(2*r^2)−(6*E*(r^2)*

SINUSOIDALY DRIVEN STRING

## Takes in the previous and the present profiles of the string and ## iterates to find the profile in the next time step. function ynext=propagate_driven(ynow,yprev,r,omega,A,n,dt) ## initialization, boundary condition ynext=ynow; ## takes in length(ynow) as an argument ## since one end of the string is driven sinuosidaly ## omega is the angular frequency of the driven force ## n is the time index, n*dt is t(n) ynext(length(ynow))=A*sin(omega*n*dt); ## entering the loop for i=2:length(ynow)-1 ynext(i) = 2*(1-r^2)*ynow(i)-yprev(i)+r^2*(ynow(i+1)+ynow(i-1)); endfor endfunction

PRIME OR NOT

## A funtion that takes a number ## in as an argument and returns ## the output p=1 if it is prime ## otherwise returns p=0. function prime_or_not(num) p=1; ## assume the number is prime. div=2; ## the smallest divisor. ## since 2*num/2> num. Forbidden Nsteps=num/2 ; ## Interval of division: the half ## of the num end itself for i= div:Nsteps ## to divide al the numbers less than ## num/2. if(rem(num,i)==0) ## any prime num cannot be ## divided with an integer. p=0; ## otherwise exit the loop. break; endif endfor if(p == 1) ## prints the result p=1 else p=0 endif endfunction

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)+(

JACOBI AND GAUSS-SIEDEL SCHEMES

##The formulas and the iteration procedure ##below are taken from the Devies Book ##Eqn 5.89 ##Constants and Initializations dx=1E-3; ##Increment in x. c1=1-2.5*dx; ## c1,c2,c3,c4 are constants in Eqn.5.89 c2=1+2.5*dx; c3=-10*dx*dx; c4=1/(2+c3); ## Set up the initial profile x=0.0:0.01:1.0; yold=0:100; hold on plot(x,yold); xlabel('x'); ylabel('y'); ##Impose the boundary conditions Nsteps=1000; yoldJac=yold; ## old y values for Jacobi Iteration Scheme yoldGauSied=yold; ## old y values for Gauss-Siedel Scheme ynewJac=yoldJac; ## new y values for Jacobi Iteration Scheme ynewGauSied=yoldGauSied; ## new y values for Gauss-Siedel Scheme for n=1:Nsteps for i=2:100 ynewJac(i)=c4*(c1*yoldJac(i+1)+c2*yoldJac(i-1)+c3*x(i)); ynewGauSied(i)=c4*(c1*yoldGauSied(i+1)+c2*ynewGauSied(i-1)+c3*x(i)); endfor yoldJac=ynewJac; yoldGauSied=ynewGauSied; endfor plot(x,ynewJac,'g*',x,ynewGauSied,'m@'); title('y vs x'); xlabel('x');

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

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

MINIMUM OF GAMMA FUNCTION

## This string finds the minimum of the Gamma function ## integral(x^(alpha-1)*exp(-x))dx) (x=[0,infinity]) ## in the interval 0<alpha<4. dx=1E-2; ## Increment in x x=0:dx:100; ## x array Gamma=[]; ## Empty Gamma function alpha=0.00:dx:4.00; ## The independent variable of Gamma function for i=1:length(alpha) ## Call the trapezoid function to ## calculate all entries of Gamma function ## corresponding to entries of alpha f=x.^(alpha(i)-1).*exp(-x); ## integrand of Gamma function Gamma(i)=trapezoid(x,f); endfor plot(alpha,Gamma); title('Gamma values vs alpha'); legend('Gamma(alpha)'); xlabel('alpha'); ylabel('Gamma'); ## The following 'for loop' finds the minimum value of Gamma ## and the corresponding alpha value. minimum=Gamma(1); for n=1:length(Gamma) if(minimum>Gamma(n)) minimum=Gamma(n); n; ## The array index where the minimum of Gamma is. m=alpha(n); ## The alpha value corresponding to the min of Gamma endi

NEWTON-RAPSON METHOD FOR HEAT FLOW

##Constants and initializations a=5.67E-8; ## Stefan-Boltzman constant[Watt/meter^2Kelvin^4] e=0.8; ## Rod surface emissivity [Dimensionless] h=20; ## Heat transfer coefficient of air flow [W/m^2-K] Tinf=Ts=25; ## Temperature of air and the walls of the close[Celcius] D=0.1; ## Diameter of the rod[meter] I2R=100; ## Electric power dissipated in rod (Ohmic Heat)[W] T=[]; ## Temperature of the rod[*C] T(1)=25; ## Initial guess of the temperature of the rod[*C] Q=[]; ## Heat function [W] Qp=[]; ## First derivative of Q wrt T [W/C*]. for i=1:100 Q(i)=pi*D*(h*(T(i)-Tinf)+e*a*(T(i)^4-Ts^4))-I2R; Qp(i)=pi*D*(h+4*e*a*T(i)^3); T(i+1)=T(i)-Q(i)/Qp(i); ## Newton-Rapson Method endfor printf('The steady state temperature is %f\n',T(i+1)) save -text HeatFlowTemp.dat ## The plot t=1:100; ##temperature for n=1:100 H(n)=pi*D*(h*(t(n)-Tinf)+e*a*(t(n)^4-Ts^4))-I2R; endfor plot(t,H) xlabel('T(Celcius)'); ylabel('Q(Watt)'); legend('Q(T)'); title(&

FACTORIAL FOR KEPLER’S MOTION

##Constants and initializations ecc=0.1; ## Eccentricity M=24.851090; ## Mean anomaly(degrees) Mr=M*pi/180; ## Mean anomaly (rads) ser=[]; ## Series in E jser=[]; ## Bessel function for m=1:100 for n=0:100 jser(n+1)=((-1)^n)*((m*ecc/2)^(2*n+m))/(factorial(n)*factorial(m+n)); jser(n+1)+=jser(n+1); endfor ser(m)=(1/m)*jser(n+1)*sin(m*Mr); ser(m)+=ser(m); endfor E=M+2*ser(m); save -text FACTORIALKEPLER.dat printf('We find E to be %f. This is the answer of Q1(c)\n',E);

BESSEL FOR KEPLER’S MOTION

##Constants and initializations ecc=0.1; ## Eccentricity M=24.851090; ## Mean anomaly(degrees) Mr=M*pi/180; ## Mean anomaly (rads) ser=[]; ## Series in E for m=1:1000 ser(m)=(1/m)*besselj(1,m*ecc)*sin(m*Mr); ser(m)+=ser(m); endfor E=M+2*ser(m); save –text BESSELKEPLER.dat printf('We find E to be %f. This is the answer of Q1(b)\n',E);

NEWTON-RAPSON METHOD FOR KEPLER’S ECCENTRICITY

##Constants and initializations ecc=0.1; ## Eccentricity M=24.851090; ## Mean anomaly(degrees) Mr=M*pi/180; ## Mean anomaly (rads) f=[]; ## The fundamental eqn we want to solve ## Kepler's Equation's zero function f(E)=E-esinE-M=0 fp=[]; ## 1st derivative of f(E). fp(E)=1+ecosE Ed=[]; ## Eccentric anomaly(degrees) Er=[]; ## Eccentric anomaly(rads) Ed(1)=M+0.85*sign(sin(Mr))*ecc; ## The eccentric anomaly initial ## guess for Danby's method(degrees) Er(1)=Ed(1)*pi/180; ## rads for i=1:100 f(i)=Ed(i)-ecc*sin(Er(i))-M; fp(i)=1+ecc*cos(Er(i)); Ed(i+1)=Ed(i)-f(i)/fp(i); ## The iteration for an updated eccentric anomaly E(i+1) ## based on a current value E(i) by Newton-Rapson Method Er(i+1)=Ed(i+1)*pi/180; endfor save -text NEWTONRAPSON.dat printf('We find E to be %f. This is the answer of Q1(a)\n',Ed(i+1));

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

NUMBER OF CLUSTER POINTS

#include <stdio.h> #include <math.h> main() { int i; double points[99][2]; FILE *fid,*write; fid=fopen("coords.txt","r"); for(i=0;i<100;i++) points[i]={x,y,z}; fscanf(fid,"%1f%1f%1f",points[i]); if(sqrt(x*x+y*y+z*z)==i) write=fopen("out.txt","w"); fprintf(fid,"%d%d\n",i,points[i]); fclose(write); else break; fclose(fid); }

INITIAL POSITIONS OF DIFFUSIVE PARTICLES

## A function takes in the number of diffusive particle ## as an argument and gets the initial positions in 2D ## as an outcome by inserting them in a square. function p=initial_positions7(Nwalkers) dt=3; ## to extend the interval of axes. a=sqrt(Nwalkers); ## step size. n=(a-1)/2; ## max range in the x and y axes. if(rem(a,2)~=1) ## Nwalkers must be odd squared. printf("The number you have entered is not an odd squared. Exiting...\n"); return endif n=(a-1)/2; ## max range in the x and y axes. x=[-n:n]'; ## integer interval in the x axis. p=zeros(a,2); ## initialize the positions of 'a' walkers. p(:,1)=x; ## equate the column of p to x vector. p=repmat(p,a,1); ## replicate positions of 'a' walkers to get a*a walkers. for i=1:a:Nwalkers ## the loop through the number walkers. m=(i+a-1)/a; ## index get the entries of x resp. p(i:i+a-1,2)=x(m) ; ## y components are a-fold degenerate. endfor plot(p(:,1),p(:,2),'b*;;') axis(dt*[-n,

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

Simple 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 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); ##Euler method n++; endwhile ##exit from the 1st while loop ##Modified Euler Method ##Constant and initializations ymid=[]; ## empty vector function evaluated at x midpoint xmid=[]; ## empty vector func. of midpoints in x ymid(1)=1; ## inital value for ymid. derymid=[]; ## derivative of y at midpoints ##Enter the 2nd while loop n=1; while (x(n)<=2) xmid(n)=x(n)+h/2; derymid(n)=xmid(n)*xmid(n)-2*ymid(

DIFFUSION-RANDOM WALK

## the function that takes in the initial positions ## of random walkers and the number of steps as ## an argument and returns the final positions ## in 2D. The figure of the walkers after the comple- ## tion of the walk also plotted. function f=diffusion(Nsteps) f=[]; Nwalkers=9; for i=1:Nsteps ## The loop through the number ## of steps each in walk. p=initial_positions(Nwalkers); ## Initial squared location of all of the ## walkers at all steps. for m=1:Nwalkers ## The loop through the desired ## number of walkers. r=rand_disc_rev(N); ## Calling the rw generator ## function which will give column vector f=[f;p+r] ## of each elements are either 1 or -1. ## Total displacement whose ## elements stems from the vector r. endfor f endfor plot(f(:,1),f(:,2),'r*;Dif;') endfunction

ISING COUPLING

## interaction energy for ising system between the spins seperated with distance r. ## J0 and alpha can be varied with respect to the ising system function J=coupling(r,J0,alpha) J=J0*exp((-r-1)/alpha); Endfunction

Equation of Motion of Proton (Coulomb Potential)

## A function of the solution of the path ## equation of a proton under the inverse ## square atraction field of an electron ## that takes in the initial seperation distance ## or the position r0 wrt center of the electron ## as an argument and returns the total time ## ttotal it takes for it to reach within 1.0m ## of the electron and plots the velocity v ## vs time graph. usage ttotal=coulomb(r0). function coulomb(r0) ## constants and initializations k=9e+9; ## coulomb field constant [Nm^2/C] q=1.6e-19; ## electronic charge [C] mp=1.7e-27; ## proton mass [m] dt=1e-4;; ## increment in time [sec] r=r0; ## initial seperation[m] t=0; ## initial time [s] v0=0; ## initial velocity [m/s] v=v0; ## 1st entry of v array [m/s] n=1; ## initialization of loop index ## since we have already n=0 in argument r0. while(r(n)>0.1); dr=v(n)*dt; ## increment that is decrease in r ## since v(n) will be negative below. r=[r;r(n)+dr]; ## decreases r in each step and ## accumulate

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