Ümit Alkuş

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((...

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 '

gnuplot> set xrange [-180:180] gnuplot> set yrange [-180:180] gnuplot>   splot asin(sin(x*pi/180)*(y-400)/ 400)*180/pi 

## 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...

## 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

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

## 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...

## 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);...

## 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...

/*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=-...

## 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

## 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...

## 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...

## 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...

## 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)*...

## 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

## 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

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

##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');...

##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 ...