Skip to main content

Posts

Showing posts from 2020

Iterative method for a nonlinear equation

Introduction: Program: ! Iterative method for a nonlinear equation real::x1,x2,check x1=-1.5 check=x1 20 check=check+1.0 x1=check 19 x2=f(x1) x1=x2 !write(*,*)check if(abs(f(x2)-x2)<=0.001)then write(*,*)"Initial choice",check+1.0 write(*,*)"The solution is" write(*,*)x1 else goto 19 endif if(x1<10.0)goto 20 end function f(q) real::q f=(q**2+6.0)/5.0 return end function Output: Summary:

gauss seidel

!PROGRAM gauss seidel integer :: i,j,k,n,no_of_ite real :: a(5,5),b(5),x(5),y(5),s,ep open(15,file="input.dat") open(16,file="output10") read(15,*)n read(15,*)no_of_ite read(15,*)ep do i=1,n read(15,*)(a(i,j),j=1,n) enddo read(15,*)(b(j),j=1,n) write(16,*)"a=" do i=1,n write(16,*)(a(i,j),j=1,n) enddo write(16,*)"b=" write(16,*)(b(j),j=1,n) do i=1,n x(i)=0.0 y(i)=0.0 enddo do k=1,no_of_ite do i=1,n s=0.0 do j=1,n if(i/=j)then s=s+a(i,j)*x(j) endif enddo x(i)=(b(i)-s)/a(i,i) enddo if((abs(x(i))-abs(y(i)))<ep)exit do i=1,n y(i)=x(i) enddo write(16,*)"no of ite=",k," x= ",(x(j),j=1,n) write(16,*)" " enddo write(16,*)"x=" write(16,*)(x(j),j=1,n) end

gauss Jocobi

!PROGRAM gauss Jocobi integer :: i,j,k,n,no_of_ite real :: a(5,5),b(5),x(5),y(5),s,ep open(15,file="input.dat") open(16,file="output10") read(15,*)n read(15,*)no_of_ite read(15,*)ep do i=1,n read(15,*)(a(i,j),j=1,n) enddo read(15,*)(b(j),j=1,n) write(16,*)"a=" do i=1,n write(16,*)(a(i,j),j=1,n) enddo write(16,*)"b=" write(16,*)(b(j),j=1,n) do i=1,n x(i)=0.0 y(i)=0.0 enddo do k=1,no_of_ite do i=1,n s=0.0 do j=1,n if(i/=j)then s=s+a(i,j)*y(j) endif enddo x(i)=(b(i)-s)/a(i,i) enddo if((abs(x(i))-abs(y(i)))<ep)exit do i=1,n y(i)=x(i) enddo write(16,*)"no of ite=",k," x= ",(x(j),j=1,n) write(16,*)" " enddo write(16,*)"x=" write(16,*)(x(j),j=1,n) end

Tridiagonal

PROGRAM Tridiagonal !The number of variable or number of equations is n. !The coefficient matrix=a(ij), i, j=1,2,3,...,n. !The lower triangular matrix=l(ij), i, j=1,2,3,...,n. !The upper matrix=u(ij), i, j=1,2,3,...,n. !b(i), i=1,2,3,...,n, the constants of the of the equations. !z(i), i=1,2,3,...,n. !x(i), i=1,2,3,...,n. integer::n,i,j,k,m real:: a(5),b(5),c(5),p(5),q(5),r(5),x(5),z(5),d(5) open(9,file="input.dat") open(10,file="output10") read(9,*)n read(9,*)(a(j),j=1,n) read(9,*)(b(j),j=1,n-1) read(9,*)(c(j),j=2,n) read(9,*)(d(j),j=1,n) write(10,*)"a="," ",(a(j),j=1,n) write(10,*)" " write(10,*)"b="," ",(b(j),j=1,n-1) write(10,*)" " write(10,*)"c="," ",(c(j),j=2,n) write(10,*)" " write(10,*)"d="," ",(d(j),j=1,n) r(1)=a(1) do i=2,n q(i)=c(i) p(i-1)=b(i-1)/r(i-1) r(i)=a(i)-p(i-1)*q(i) write(10,*)"i=",i,"p,q,r=",p(i),q(i),r(i) enddo writ

LU-Decomposition

! PROGRAM LU-Decomposition !The number of variable or number of equations is n. !The coefficient matrix=a(ij), i, j=1,2,3,...,n. !The lower triangular matrix=l(ij), i, j=1,2,3,...,n. !The upper matrix=u(ij), i, j=1,2,3,...,n. !b(i), i=1,2,3,...,n, the constants of the of the equations. !z(i), i=1,2,3,...,n. !x(i), i=1,2,3,...,n. integer::n,i,j,k,m real:: a(5,5),l(5,5),u(5,5),b(5),x(5),z(5),ratio,s open(3,file="input.dat") open(4,file="output2") read(3,*)n do i=1,n read(3,*)(a(i,j),j=1,n) enddo write(4,*) "A=" do i=1,n write(4,*)(a(i,j),j=1,n) enddo write(4,*)" " read(3,*)(b(j),j=1,n) write(4,*)"B=",(b(j),j=1,n) write(4,*)" " do i=1,n l(i,1)=a(i,1) u(i,i)=1.0 enddo do j=2,n u(1,j)=a(1,j)/l(1,1) enddo do k=2,n !To determine columns of l do i=k,n s=0 do j=1,i-1 s=s+l(i,j)*u(j,k) enddo l(i,k)=a(i,k)-s enddo !!To determine rows of u do i=k+1,n if(i/=k)then s=0 do j=1,i-1 s=s+l(k,j)*u(j,i) enddo endif u(k,i)=(a(k,i)-s)/l(k,k) enddo

Gauss Jordan

PROGRAM Gauss_Jordan ! Gauss Jordan method to find the solution of the system of linear equations ! The number of equations or the number of variables of the system ! of linear equations is n ! The coefficient matrix =a(ij) : i, j=1,2,3,...,n ! The constant vector =b(i) : i=1,2,3,...,n ! The vector of the variables =x(i) : i=1,2,3,...,n integer::n,i,j,k,m real:: a(5,5),b(5),x(5),r open(3,file="input.dat") open(4,file="output2") read(3,*)n do i=1,n read(3,*)(a(i,j),j=1,n) enddo read(3,*)(b(j),j=1,n) write(4,*) "a=" do i=1,n write(4,*)(a(i,j),j=1,n) enddo write(4,*)" " write(4,*)"b=",(b(j),j=1,n) do k=1,n do i=1,n if(i/=k)then r=a(i,k)/a(k,k) do j=k,n a(i,j+1)=a(i,j+1)-a(k,j+1)*r enddo b(i)=b(i)-b(k)*r endif enddo enddo do i=1,n x(i)=b(i)/a(i,i) enddo write(4,*)" " write(4,*)"The solution set is ","",(x(j),j=1,n),"" end PROGRAM

Gauss Elimination

PROGRAM Gauss_Elimination ! n is number of variables or number of equations ! a[ij], i, j=1,2,3,...,n is coefficient matrix, b[i], i=1,2,3,...,n is constant ! vector, x[i], i=1,2,3...,n is the vector of variables integer::n,i,j,k real:: a(5,5),b(5),x(5),ratio,s open(1,file="input.dat") open(2,file="output.dat") read(1,*)n write(2,*)"The number of variables or number of equations, n=",n do i=1,n read(1,*)(a(i,j),j=1,n) enddo write(2,*)"The cefficient matrix, A=" do i=1,n write(2,*)(a(i,j),j=1,n) enddo do j=1,n read(1,*)b(j) enddo write(2,*)"The constant vector, B=" write(2,*)(b(j),j=1,n) do k=1,n-1 do i=k+1,n ratio=a(i,k)/a(k,k) do j=k,n a(i,j)=a(i,j)-ratio*a(k,j) enddo b(i)=b(i)-ratio*b(k) enddo enddo write(2,*) do i=1,n write(2,*)(a(i,j),j=1,n) enddo x(n)=b(n)/a(n,n) do i=n-1,1,-1 s=0.0 do j=n,i+1,-1 s=s+a(i,j)*x(j) enddo x(i)=(b(i)-s)/a(i,i) enddo write(2,*) write(2,*)(x(j),j=1,n) end PROGRAM

Damped harmonic oscillator using Runge Kutta method

PROGRAM Damped_harmonic_oscillator integer::i,n,k real::t,dt,a,x,v,k1,k2,k3,k4 f(t,v)=-a*v-x f1(t,x)=v open(15,FILE="p2out3.dat",STATUS="unknown") a=1.0 t=0;x=1;v=1;n=300 dt=0.1 write(15,*)t,x do j=1,n k1=dt*f(t,v) k2=dt*f(t+0.5*dt,v+0.5*k1) k3=dt*f(t+0.5*dt,v+0.5*k2) k4=dt*f(t+dt,v+k3) v=v+k1/6+k2/3+k3/3+k4/6 k1=dt*f1(t,x) k2=dt*f1(t+0.5*dt,x+0.5*k1) k3=dt*f1(t+0.5*dt,x+0.5*k2) k4=dt*f1(t+dt,x+k3) x=x+k1/6+k2/3+k3/3+k4/6 t=t+dt write(15,*)t,x enddo END PROGRAM

Solve second order ODE of One-dimensional harmonic oscillator using Euler Method

!PROGRAM to solve second order ODE of the form d^2xdt^2 + x = 0 !One-dimensional harmonic oscillator using Euler Method real::t,v,x,dt integer::time=20 open(1,file="harmonic_oscillator.out") open(2,file="harmonic_oscillator1.out") t=0;x=0.1;v=0;dt=0.1 write(1,*)t,x 1 v=v-dt*x x=x+dt*v if(t<=time) then t=t+dt write(1,*)t,x write(2,*)x,v go to 1 endif end

Predictor-Corrector method

real::lambda,t,N,dt,NP,NC,N1,t1 f(t,N)=-lambda*N open(1,file="output.dat") t=0;N=0.5;lambda=0.1;dt=0.1 write(1,*)t,N N1=N+dt*f(t,N) t1=t+dt write(1,*)t1,N1 2 NP=N1+0.5*dt*(f(t,N)+f(t1,N1)) NC=N1+0.5*dt*(f(t,N)+f(t1,NP)) if(t1<=100) then N=N1 N1=NP t=t1 t1=t1+dt write(1,*)t1,NC goto 2 endif end

Range-Kutta method to solve di/dt+iR/L=V/L

! PROGRAM RK method to solve di/dt+iR/L=V/L real::t,i,k1,k2,k3,k4,V,L,R,io integer::j,jj,n f(t,i)=V/L-R/L*i V=20; L=0.4; R=2 open(1,file="RK_method.out") open(2,file="RK_analytical.out") t=0;i=0;h=0.4;n=20 write(1,*)t,i do j=1,n k1=h*f(t,i) k2=h*f(t+0.5*h,i+0.5*k1) k3=h*f(t+0.5*h,i+0.5*k2) k4=h*f(t+h,i+k3) i=i+k1/6+k2/3+k3/3+k4/6 t=t+h write(1,*)t,i enddo io=V/R t=0 do jj=1,n i=io*(1-exp(-R/L*t)) write(2,*)t,i t=t+h enddo end

Modified Euler method to solve m*dv/dt = mg - bv with v(0)=0

! PROGRAM Modified Euler method to solve m dv/dt= mg- bv with v(0)=0 real::t,v,g,b,m,t1,v1 f(t,v)=g-(b/m)*v open(3,file="Mod_euler.out") m=1.0 g=9.8 b=10.0 t=0;v=0;h=0.5 1 write(3,*)t,v v1=v+h*f(t,v) t1=t+h v=v+0.5*h*(f(t,v)+f(t1,v1)) t=t1 if(t<=5) goto 1 open(4,file="analytical.out") t=0;v=0 2 write(4,*)t,v t=t+h v=(m*g)/b*(1-exp(-(b/m)*t)) if(t<=5)goto 2 end

Non-linear: Lorentz attractor using Euler’s method

PROGRAM Lorentz_attractor PARAMETER(kx=100000) INTEGER::i,j,n,k REAL::a,r,b,t(kx),dt,x(kx),y(kx),z(kx) !a=sigma OPEN(15,FILE="p3.dat",STATUS="unknown") r=28.0 a=10.0 b=8.0/3.0 x(1)=0.0 !Initial points y(1)=1.0 z(1)=0.0 dt=(0.001) t(1)=0.0 DO i=1,kx-1 x(i+1)=x(i)+dt*a*(y(i)-x(i)) y(i+1)=y(i)+dt*(r*x(i)-y(i)-x(i)*z(i)) z(i+1)=z(i)+dt*(x(i)*y(i)-b*z(i)) t(i+1)=t(i)+dt ENDDO DO i=1,kx-1 WRITE(15,1)t(i),x(i),y(i),z(i) ENDDO 1 FORMAT(4(F7.3,1x)) END PROGRAM

Non-linear: Rabbit Sheep using Euler’s method

PROGRAM Rabbit_Sheep PARAMETER(kx=10000) INTEGER::i,j,n,k REAL::t(kx),dt,x(50,kx),y(50,kx) OPEN(25,FILE="p6in.dat") OPEN(15,FILE="p6out.dat",STATUS="unknown") READ(25,*) k DO n=1,k READ(25,*) x(n,1),y(n,1) ENDDO dt=0.001 t(1)=0.0 DO n=1,k DO i=1,kx-1 x(n,i+1)=x(n,i)+dt*(x(n,i)*(3.0-x(n,i)-2.0*y(n,i))) y(n,i+1)=y(n,i)+dt*(y(n,i)*(2.0-x(n,i)-y(n,i))) t(i+1)=t(i)+dt ENDDO ENDDO DO i=1,kx-1 WRITE(15,1)t(i),(x(n,i),n=1,k),(y(n,i),n=1,k) ENDDO 1 FORMAT(45(F7.3,1x)) END PROGRAM

Non-linear: Van der Pol Oscillator using Euler’s method

PROGRAM van_der_Pol_Oscillator parameter (kx=50000) integer::i,j,n,k real::t(kx),x(50,kx),v(50,kx),u open(15,file="vpolin.dat",status="unknown") open(25,file="vpolout.dat",status="unknown") write(*,*) "enter the value of mue" read(*,*) u read(15,*) k dt=0.001 t(1)=0.0 do n=1,k read(15,*) x(n,1),v(n,1) enddo do n=1,k do i=1,kx-1 x(n,i+1)=x(n,i)+dt*v(n,i) v(n,i+1)=v(n,i)-dt*u*((x(n,i)**2-4)*v(n,i))-dt*x(n,i) t(i+1)=t(i)+dt enddo enddo do i=1,kx write(25,3) t(i),(x(n,i),n=1,k),(v(n,i),n=1,k) enddo 3 format (20(f7.3,1x)) end

Pendulum motion with nonlinearlity

Program Pendulum_motion_with_nonlinearlity implicit none real::x,y,dt real,parameter::pi=3.141593 integer::i,n open(11,file="ho1.dat",status="unknown") open(12,file="ho2.dat",status="unknown") open(13,file="ho3.dat",status="unknown") open(14,file="ho4.dat",status="unknown") open(15,file="ho5.dat",status="unknown") open(16,file="ho6.dat",status="unknown") !Initial condition y=0 and x=.5 x=0.5; y=0.0; dt=0.01 write(11,*)x,y do i=1,1500 x=x + y*dt y=y - dt*sin(x) write(11,*)x,y enddo close(11) !Initial condition y=0 x=1.8 x=1.8; y=0.0; dt=0.01 write(12,*)x,y do i=1,1500 x=x + y*dt y=y - dt*sin(x) write(12,*)x,y enddo close(12) !Initial condition x=pi-0.0001, y=-0.0001 x=pi-0.0001; y=-0.0001; dt=0.01 write(13,*)x,y do i=1,1500 x=x + y*dt y=y - dt*sin(x) write(13,*)x,y enddo close(13) !Initial condition x=-pi+0.0001 y=0.0001end x=-pi+0.0001; y=0.0001; dt=0.01 write(14,*)x