MHB Gauss-Siedel Matrix to solve Elliptic Equation

[ra_forever8]

Qs: Write a FORTRAN program to approximately solve elliptic equation : - u_xx - u_yy=1 with the five-point finite difference formula in the right-angled triangle, sides 1,2,3^(1/2) using the Gauss-Seidel matrix solver withe the boundary conditions u=0 on all sides. Take a mesh spacing to be h=0.01. (You'll need linear interpolation on the long side.)

=>

  PROGRAM Gauss_Seidel
IMPLICIT NONE

! Declare Variables
Real, allocatable :: x(:),y(:),u(:,:),u_old(:,:)
Real:: h,tolerence,error
Integer:: i,j,JI,NI

h=0.01
JI=100
NI=173
error = 1.d0
tolerence = 10E-4
 ! Total number of space stepsize
allocate (x(0:JI),y(0:NI),u(0:JI+1,0:NI+1),u_old(0:JI+1,0:NI+1))
open(10,file='Gauss_Seidel.m')  !Opening files in Matlab

!Initial Conditions
x(0)= 0
x(JI)= 1.0
y(0)= 0
y(NI)= SQRT(3.0)

do i=0,JI
do j=0,NI
x(i)= i*h   ! x-axix, x starts from 0 to 1 
y(j)= j*h   ! y-axis  y starts from 0 to SQRT(3.0)
u(i,j)= 0         ! Entire Boundary is zero
end do
end do

while (error .GT. tolerence) do  ! To stop
  do i=1, JI-1
    do j=1,NI-1
u_old(i,j)= u(i,j)  ! To [U]store[/U] the old values

!Using 5-point scheme Formulae and rearranging the equation
u(i,j)= 0.25*(u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1)+h**2) 
     end do
     end do

error =0.d0        ! Now, error reading the value of zero
     do i=1,JI-1
       do j=1, NI-1
error = error + abs(u(i,j)- u_old(i,j))  ! To Stop
        end do
         end do
        end do

 do i=1, JI-1
    do j=1,NI-1

 u(i+1,j)= - (alpha * u(i,j))/ beta    ! lies outside the long side
 u(i,j+1)= - (a * u(i,j))/ b           ! lies outside the long side
 
     end do
     end do
!Print out the Approximate solution in MATLAB to get output and plot 
write(10,*)  'x=['                           
     do i=0, JI
       write(10,*) x(i)
 end do
write(10,*) ']'

write(10,*)  'y=['                           
       do j=0,NI
write(10,*) y(j)
 end do
write(10,*) ']'

write(10,*)  'u=['                           
     do i=0, JI
       do j=0,NI
write(10,*) u(i,j)
 end do
 end do
write(10,*) ']'

write(10,*) "[X,Y] = meshgrid(x,y)"              !Ploting diagram x,y,u
write(10,*) "Z=reshape(u,length(y), length(x))"
write(10,*) "contour(X,Y,Z)"
write(10,*) "xlabel('x'),ylabel('y'),legend('Approximate Gauss Seidel')"
close(10) 

END PROGRAM Gauss_Seidel

(Note: THIS IS WHAT I HAVE GOT SO FAR, STILL LITTLE BIT INCOMPLETE.
THE LENGTH OF THE LONG SIDE OF THE TRIANGLE IS 2. SINCE THE POINTS ON THE LONG SIDE DO NOT COINCIDE WITH THE GRID POINTS, THE POINTS u(i,j+1) and u(i+1,j) lies outside the long side of the triangle, the distance between the points u(i,j) to u(i,j+1) and u(i,j) to u(i+1,j) are zero which is known using the linear interpolation. Alpha is the distance from zero to point u(i + 1, j) and beta is the distance from point u(i, j) to zero.
Similarly, a is the distance from zero to point u(i, j+1) and b is the distance from point u(i, j) to zero.
The points u(i -1, j) and u(i, j - 1) lies inside the long side of the triangle, the distance from these points to u(i, j) remains unchanged and does not affect in the five-point scheme.
Can anyone help me after this. What can do with alpha, beta,a and b? Do I need to know values for alpha,beta, a and b?
Please help me )

 

[member38644]

In order to complete the program, you will need to calculate the values for alpha, beta, a, and b. These values represent the distances from the grid points to the long side of the triangle, and they will affect the calculation of the five-point scheme.

To calculate these values, you can use linear interpolation. This involves finding the distance between the grid points and the long side, and then using that distance to calculate the distances to the points u(i+1,j) and u(i,j+1).

Once you have calculated these values, you can use them in the equations for u(i+1,j) and u(i,j+1) in the program. This will ensure that the calculation takes into account the points that lie outside the long side of the triangle.

You can also use these values to plot the approximate solution in MATLAB, as shown in the program. This will help you visualize the solution and make any necessary adjustments to the program.

Overall, the key to completing this program is to accurately calculate the values for alpha, beta, a, and b. Once you have these values, you can use them in the program and plot the solution to see if it matches the expected solution for the given boundary conditions and equation.