Jacobi iterative method to solve elliptic equation
Jacobi iterative method to solve elliptic equation
(OP)
Qs: Write a FORTRAN program to approximately solve elliptic equation :-u_xx-u_yy=1 on a unit square -1<x<1, -1<y<1, with u=0 on the entire boundary.
Use the five -point finite difference formula on a uniform mesh of size h in each direction.
Set up the matrix system Au=b (taking the natural ordering of the nodes) and program the Jacobi iteration method to solve it, i.e.
u^k = (I-D^-1 A)u^(k-1) + D^-1 b, where k is the iteration index and D is the diagonal of A.
=> So far my program
PROGRAM Jacobi_method
IMPLICIT NONE
! Declare Variables
Real, allocatable :: x(:),u(:,:)
Real:: m,xy,h,Smax
Integer:: i,j,JI
Print *, 'Enter the space size:'
read*, xy
Print*, 'Enter the final space:'
read*, Smax
h=Smax/xy !The size of spacestep
Print*, 'This gives stepsize of space h=',h
JI=20 ! Total number of space stepsize
allocate (x(0:JI),u(0:JI+1,0:JI+1))
open(10,file='Jacobi_method.m')
do i=1, JI
do j=1,JI
! Using 5-point scheme Formulae
u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1)-4*u(i,j)=-h**2
end do
end do
there are alots of code missing .can anyone please help me. Iam really stuck in this programing.....Thanks alot in advance.
END PROGRAM Jacobi_method
Use the five -point finite difference formula on a uniform mesh of size h in each direction.
Set up the matrix system Au=b (taking the natural ordering of the nodes) and program the Jacobi iteration method to solve it, i.e.
u^k = (I-D^-1 A)u^(k-1) + D^-1 b, where k is the iteration index and D is the diagonal of A.
=> So far my program
PROGRAM Jacobi_method
IMPLICIT NONE
! Declare Variables
Real, allocatable :: x(:),u(:,:)
Real:: m,xy,h,Smax
Integer:: i,j,JI
Print *, 'Enter the space size:'
read*, xy
Print*, 'Enter the final space:'
read*, Smax
h=Smax/xy !The size of spacestep
Print*, 'This gives stepsize of space h=',h
JI=20 ! Total number of space stepsize
allocate (x(0:JI),u(0:JI+1,0:JI+1))
open(10,file='Jacobi_method.m')
do i=1, JI
do j=1,JI
! Using 5-point scheme Formulae
u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1)-4*u(i,j)=-h**2
end do
end do
there are alots of code missing .can anyone please help me. Iam really stuck in this programing.....Thanks alot in advance.
END PROGRAM Jacobi_method
RE: Jacobi iterative method to solve elliptic equation
CODE -->
I have added some code, but i guess there is still some code missing to answer the question
With the above code,
I have got Run-time error: Undefined variables element or arrays function in this line : u(i,j)= 0.25*(u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1)+h**2)
Please anyone help me
RE: Jacobi iterative method to solve elliptic equation
- you must solve a problem with -1 < x < +1 and -1 < y < -1 with boundary conditions and a constant step size common to both directions
I suggest to generalize immediately
xmin < x < xmax and ymin < y < ymax with two different step sizes hx and hy
if nx is the number of intervals in x direction, then you must construct and array X with nx+1 values, for instance X(0:nx) verifying :
X(0)==xmin ; x(nx) == xmax => X(i+1)=xmin+(xmax-xmin)/nx*i for i from 0 to nx
of course hx=(xmax-xmin)/nx.
I suggest that the main input data should be xmin,xmax,ymin,ymax,nx and ny
- your partial derivatives becomes :
u_xx(i,j)=(u(i+1,j)-2*u(i,j)+u(i-1,j))/hx**2
u_yy(i,j)=(u(i,j+1)-2*u(i,j)+u(i,j-1))/hy**2
please, notice that I use different notations between u and u_xx and u_yy. In your programming, there is a full confusion between
the laplacian L=u_xx+u_yy and the function u. In addition, the problem to solve is L=1 (I suggest again to generalize to L=val, val being
a real value read from the standard input like xmin, xmax, ymin xmax, nx and ny.
I have the feeling that you try to solve u=L. This is not the right problem !
- How many unknows do you have ? Clearly, you have not the right answer !
the unknows are u(i,j) from i=1 to i=nx-1 and j=1 to j=ny-1 => (nx-1)*(ny-1) values
to manage boundary conditions easily, the array u may be declared u(0:nx,0:ny). But the derivatives u_xx can be computed only for i from 1 to nx-1 and
u_yy from j from 1 to ny-1 => L=u_xx+u_yy has sense only for 1<=i<=nx-1 and 1<=j <=ny-1
=> Please : deduce the right declaration for the array L representing u_xx+u_yy
- The text tells a story about an iterative algorithm, in particular a version u^k and another u^(k-1). How many versions of arrays u do you need ?
- The text mention a diagonal. How do you declare it, How do you construct it ?
François Jacq
RE: Jacobi iterative method to solve elliptic equation
Qs: Write a FORTRAN program to approximately solve elliptic equation :-u_xx-u_yy=1 on a unit square -1<x<1, -1<y<1, with u=0 on the entire boundary.
Use the five -point finite difference formula on a uniform mesh of size h in each direction.
Set up the matrix system Au=b (taking the natural ordering of the nodes) and program the Jacobi iteration method to solve it, i.e.
u^k = (I-D^-1 A)u^(k-1) + D^-1 b, where k is the iteration index and D is the diagonal of A.
=>
What iam trying to solve is The poisson euqation :-u_xx-u_yy= 1 not laplace equations, by using Jacobi iterative method for solving linear systems such as
A*u= b
u= A^(-1)*b
Where b is the right hand side, u is unknown and A is matrix
For this, we use a sequence u^(k) which converges to fixed point (solution) u.
u^k = (I-D^-1 A)u^(k-1) + D^-1 *b
where I is the identity matrix ,k is the iteration index and D is the diagonal of A.
In the jacobi's method, we choose M=D, A=M-N, and N= E+F where M is an invertible matrix.
u^(k+1)= D^-1(E+F)u^(k)+D^-1*b
Where E is the strcitly lower triangular part of A
F is the strcitly upper triangular part of A
I have to solve this problem Using 5-point scheme Formulae: u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1)-4*u(i,j)=-h^2
And For your question about different size of h for x and y direction. I think h will be same for both directions beacause -1<x<1, -1<y<1 and uniform mesh size of h is equal in each direction.
And also, your partial derivatives :
u_xx(i,j)=(u(i+1,j)-2*u(i,j)+u(i-1,j))/hx**2
u_yy(i,j)=(u(i,j+1)-2*u(i,j)+u(i,j-1))/hy**2
Using the 5-point scheme Formulae, I have derived your partial derivatives and since h is same both x and y
i got this equations u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1)-4*u(i,j)=-h^2
I hope i have make it clear, if not please reply me back.
This is what iam trying to solve but i know my fortran code are not correct and some of the coding are missing too. Please help me if you understand the problem.
RE: Jacobi iterative method to solve elliptic equation
And same for the poisson euqation : - u_xx - u_yy = 1
RE: Jacobi iterative method to solve elliptic equation
! The Entire Boundary Conditions
u(0,1)=0 !West u(i-1,j), x direction
u(2,1)=0 ! East u(i+1,j), x direction
u(1,2)=0 !North u(i,j+1), y direction
u(1,0)=0 !South u(i,j-1), y direction
because The Entire Boundary Conditions of u=0, when i=1, j=1
u(0,1) ,West u(i-1,j), x direction
u(2,1), East u(i+1,j), x direction
u(1,2), North u(i,j+1), y direction
u(1,0), South u(i,j-1), y direction
are unknown and are not equal to zero because they are single point of grid inside the unit square.
West u(i-1,j)and East u(i+1,j): x direction, North u(i,j+1) and South u(i,j-1):y direction. Between East,West,North and South there is a point u(i,j) in the middle which is unknown and thats what i need to find out.
I wish i could draw the diagram to clarify easy.
RE: Jacobi iterative method to solve elliptic equation
I have read your problem carefully. But I don't have to solve it : this is YOUR TASK !
For piece of information, the operator L transforming u into u_xx+u_yy is called the Laplace operator. I never spoke about the Laplace equation which is L(u)==0 and is a particular case of the Poisson's equation.
Rather than sticking to the exact text of a problem, especially in the programming area, is often much more clear and powerful to generalize a little bit. For instance, I proposed to solve L=val : choosing val=-1. was corresponding to your exact problem.
You did not answer most of my questions... If you answer them correctly, you should be able to build your program easily without any additional help.
François Jacq
RE: Jacobi iterative method to solve elliptic equation
Refering back to qs agian:
Qs: Write a FORTRAN program to approximately solve elliptic equation : - u_xx - u_yy= 1 on a unit square -1<x<1, -1<y<1, with u=0 on the entire boundary.
Use the five -point finite difference formula on a uniform mesh of size h in each direction.
Set up the matrix system Au=b (taking the natural ordering of the nodes) and program the Jacobi iteration method to solve it, i.e.
u^k = (I-D^-1 A)u^(k-1) + D^-1 b, where k is the iteration index and D is the diagonal of A.
=> My new code
CODE
There are no errors and i manage to generate some outputs.
RE: Jacobi iterative method to solve elliptic equation
ApproximateSolution =
x-axis y-axis u(i,j)
[ -1.00000 -1.00000 0.00000
-0.900000 -0.900000 2.500000E-03
-0.800000 -0.900000 3.125000E-03
-0.700000 -0.900000 3.281250E-03
-0.600000 -0.900000 3.320313E-03
-0.500000 -0.900000 3.330078E-03
-0.400000 -0.900000 3.332520E-03
-0.300000 -0.900000 3.333130E-03
-0.200000 -0.900000 3.333283E-03
-9.999999E-02 -0.900000 3.333321E-03
1.490116E-08 -0.900000 3.333330E-03
0.100000 -0.900000 3.333333E-03
0.200000 -0.900000 3.333333E-03
0.300000 -0.900000 3.333333E-03
0.400000 -0.900000 3.333333E-03
0.500000 -0.900000 3.333333E-03
0.600000 -0.900000 3.333333E-03
0.700000 -0.900000 3.333333E-03
0.800000 -0.900000 3.333333E-03
0.900000 -0.900000 3.333333E-03
-0.900000 -0.800000 3.125000E-03
-0.800000 -0.800000 4.062500E-03
-0.700000 -0.800000 4.335938E-03
-0.600000 -0.800000 4.414062E-03
-0.500000 -0.800000 4.436035E-03
-0.400000 -0.800000 4.442139E-03
-0.300000 -0.800000 4.443817E-03
-0.200000 -0.800000 4.444275E-03
-9.999999E-02 -0.800000 4.444399E-03
1.490116E-08 -0.800000 4.444432E-03
0.100000 -0.800000 4.444441E-03
0.200000 -0.800000 4.444444E-03
0.300000 -0.800000 4.444445E-03
0.400000 -0.800000 4.444445E-03
0.500000 -0.800000 4.444445E-03
0.600000 -0.800000 4.444445E-03
0.700000 -0.800000 4.444445E-03
0.800000 -0.800000 4.444445E-03
0.900000 -0.800000 4.444445E-03
-0.900000 -0.700000 3.281250E-03
-0.800000 -0.700000 4.335938E-03
-0.700000 -0.700000 4.667969E-03
-0.600000 -0.700000 4.770508E-03
-0.500000 -0.700000 4.801636E-03
-0.400000 -0.700000 4.810944E-03
-0.300000 -0.700000 4.813690E-03
-0.200000 -0.700000 4.814492E-03
-9.999999E-02 -0.700000 4.814723E-03
1.490116E-08 -0.700000 4.814789E-03
0.100000 -0.700000 4.814808E-03
0.200000 -0.700000 4.814813E-03
0.300000 -0.700000 4.814814E-03
0.400000 -0.700000 4.814815E-03
0.500000 -0.700000 4.814815E-03
0.600000 -0.700000 4.814815E-03
0.700000 -0.700000 4.814815E-03
0.800000 -0.700000 4.814815E-03
0.900000 -0.700000 4.814815E-03
-0.900000 -0.600000 3.320313E-03
-0.800000 -0.600000 4.414062E-03
-0.700000 -0.600000 4.770508E-03
-0.600000 -0.600000 4.885254E-03
-0.500000 -0.600000 4.921723E-03
-0.400000 -0.600000 4.933167E-03
-0.300000 -0.600000 4.936714E-03
-0.200000 -0.600000 4.937802E-03
-9.999999E-02 -0.600000 4.938131E-03
1.490116E-08 -0.600000 4.938230E-03
0.100000 -0.600000 4.938260E-03
0.200000 -0.600000 4.938268E-03
0.300000 -0.600000 4.938271E-03
0.400000 -0.600000 4.938271E-03
0.500000 -0.600000 4.938272E-03
0.600000 -0.600000 4.938272E-03
0.700000 -0.600000 4.938272E-03
0.800000 -0.600000 4.938272E-03
0.900000 -0.600000 4.938272E-03
-0.900000 -0.500000 3.330078E-03
-0.800000 -0.500000 4.436035E-03
-0.700000 -0.500000 4.801636E-03
-0.600000 -0.500000 4.921723E-03
-0.500000 -0.500000 4.960862E-03
-0.400000 -0.500000 4.973507E-03
-0.300000 -0.500000 4.977555E-03
-0.200000 -0.500000 4.978839E-03
-9.999999E-02 -0.500000 4.979243E-03
1.490116E-08 -0.500000 4.979368E-03
0.100000 -0.500000 4.979407E-03
0.200000 -0.500000 4.979419E-03
0.300000 -0.500000 4.979423E-03
0.400000 -0.500000 4.979424E-03
0.500000 -0.500000 4.979424E-03
0.600000 -0.500000 4.979424E-03
0.700000 -0.500000 4.979424E-03
0.800000 -0.500000 4.979424E-03
0.900000 -0.500000 4.979424E-03
-0.900000 -0.400000 3.332520E-03
-0.800000 -0.400000 4.442139E-03
-0.700000 -0.400000 4.810944E-03
-0.600000 -0.400000 4.933167E-03
-0.500000 -0.400000 4.973507E-03
-0.400000 -0.400000 4.986754E-03
-0.300000 -0.400000 4.991077E-03
-0.200000 -0.400000 4.992479E-03
-9.999999E-02 -0.400000 4.992931E-03
1.490116E-08 -0.400000 4.993075E-03
0.100000 -0.400000 4.993121E-03
0.200000 -0.400000 4.993135E-03
0.300000 -0.400000 4.993140E-03
0.400000 -0.400000 4.993141E-03
0.500000 -0.400000 4.993142E-03
0.600000 -0.400000 4.993142E-03
0.700000 -0.400000 4.993142E-03
0.800000 -0.400000 4.993142E-03
0.900000 -0.400000 4.993142E-03
-0.900000 -0.300000 3.333130E-03
-0.800000 -0.300000 4.443817E-03
-0.700000 -0.300000 4.813690E-03
-0.600000 -0.300000 4.936714E-03
-0.500000 -0.300000 4.977555E-03
-0.400000 -0.300000 4.991077E-03
-0.300000 -0.300000 4.995539E-03
-0.200000 -0.300000 4.997005E-03
-9.999999E-02 -0.300000 4.997484E-03
1.490116E-08 -0.300000 4.997640E-03
0.100000 -0.300000 4.997690E-03
0.200000 -0.300000 4.997707E-03
0.300000 -0.300000 4.997712E-03
0.400000 -0.300000 4.997713E-03
0.500000 -0.300000 4.997714E-03
0.600000 -0.300000 4.997714E-03
0.700000 -0.300000 4.997714E-03
0.800000 -0.300000 4.997714E-03
0.900000 -0.300000 4.997714E-03
-0.900000 -0.200000 3.333283E-03
-0.800000 -0.200000 4.444275E-03
-0.700000 -0.200000 4.814492E-03
-0.600000 -0.200000 4.937802E-03
-0.500000 -0.200000 4.978839E-03
-0.400000 -0.200000 4.992479E-03
-0.300000 -0.200000 4.997005E-03
-0.200000 -0.200000 4.998502E-03
-9.999999E-02 -0.200000 4.998997E-03
1.490116E-08 -0.200000 4.999159E-03
0.100000 -0.200000 4.999212E-03
0.200000 -0.200000 4.999230E-03
0.300000 -0.200000 4.999235E-03
0.400000 -0.200000 4.999237E-03
0.500000 -0.200000 4.999238E-03
0.600000 -0.200000 4.999238E-03
0.700000 -0.200000 4.999238E-03
0.800000 -0.200000 4.999238E-03
0.900000 -0.200000 4.999238E-03
-0.900000 -9.999999E-02 3.333321E-03
-0.800000 -9.999999E-02 4.444399E-03
-0.700000 -9.999999E-02 4.814723E-03
-0.600000 -9.999999E-02 4.938131E-03
-0.500000 -9.999999E-02 4.979243E-03
-0.400000 -9.999999E-02 4.992931E-03
-0.300000 -9.999999E-02 4.997484E-03
-0.200000 -9.999999E-02 4.998997E-03
-9.999999E-02 -9.999999E-02 4.999498E-03
1.490116E-08 -9.999999E-02 4.999665E-03
0.100000 -9.999999E-02 4.999720E-03
0.200000 -9.999999E-02 4.999737E-03
0.300000 -9.999999E-02 4.999743E-03
0.400000 -9.999999E-02 4.999745E-03
0.500000 -9.999999E-02 4.999746E-03
0.600000 -9.999999E-02 4.999746E-03
0.700000 -9.999999E-02 4.999746E-03
0.800000 -9.999999E-02 4.999746E-03
0.900000 -9.999999E-02 4.999746E-03
-0.900000 1.490116E-08 3.333330E-03
-0.800000 1.490116E-08 4.444432E-03
-0.700000 1.490116E-08 4.814789E-03
-0.600000 1.490116E-08 4.938230E-03
-0.500000 1.490116E-08 4.979368E-03
-0.400000 1.490116E-08 4.993075E-03
-0.300000 1.490116E-08 4.997640E-03
-0.200000 1.490116E-08 4.999159E-03
-9.999999E-02 1.490116E-08 4.999665E-03
1.490116E-08 1.490116E-08 4.999832E-03
0.100000 1.490116E-08 4.999888E-03
0.200000 1.490116E-08 4.999906E-03
0.300000 1.490116E-08 4.999912E-03
0.400000 1.490116E-08 4.999915E-03
0.500000 1.490116E-08 4.999915E-03
0.600000 1.490116E-08 4.999916E-03
0.700000 1.490116E-08 4.999916E-03
0.800000 1.490116E-08 4.999916E-03
0.900000 1.490116E-08 4.999916E-03
-0.900000 0.100000 3.333333E-03
-0.800000 0.100000 4.444441E-03
-0.700000 0.100000 4.814808E-03
-0.600000 0.100000 4.938260E-03
-0.500000 0.100000 4.979407E-03
-0.400000 0.100000 4.993121E-03
-0.300000 0.100000 4.997690E-03
-0.200000 0.100000 4.999212E-03
-9.999999E-02 0.100000 4.999720E-03
1.490116E-08 0.100000 4.999888E-03
0.100000 0.100000 4.999944E-03
0.200000 0.100000 4.999963E-03
0.300000 0.100000 4.999969E-03
0.400000 0.100000 4.999971E-03
0.500000 0.100000 4.999971E-03
0.600000 0.100000 4.999972E-03
0.700000 0.100000 4.999972E-03
0.800000 0.100000 4.999972E-03
0.900000 0.100000 4.999972E-03
-0.900000 0.200000 3.333333E-03
-0.800000 0.200000 4.444444E-03
-0.700000 0.200000 4.814813E-03
-0.600000 0.200000 4.938268E-03
-0.500000 0.200000 4.979419E-03
-0.400000 0.200000 4.993135E-03
-0.300000 0.200000 4.997707E-03
-0.200000 0.200000 4.999230E-03
-9.999999E-02 0.200000 4.999737E-03
1.490116E-08 0.200000 4.999906E-03
0.100000 0.200000 4.999963E-03
0.200000 0.200000 4.999981E-03
0.300000 0.200000 4.999988E-03
0.400000 0.200000 4.999990E-03
0.500000 0.200000 4.999991E-03
0.600000 0.200000 4.999991E-03
0.700000 0.200000 4.999991E-03
0.800000 0.200000 4.999991E-03
0.900000 0.200000 4.999991E-03
-0.900000 0.300000 3.333333E-03
-0.800000 0.300000 4.444445E-03
-0.700000 0.300000 4.814814E-03
-0.600000 0.300000 4.938271E-03
-0.500000 0.300000 4.979423E-03
-0.400000 0.300000 4.993140E-03
-0.300000 0.300000 4.997712E-03
-0.200000 0.300000 4.999235E-03
-9.999999E-02 0.300000 4.999743E-03
1.490116E-08 0.300000 4.999912E-03
0.100000 0.300000 4.999969E-03
0.200000 0.300000 4.999988E-03
0.300000 0.300000 4.999994E-03
0.400000 0.300000 4.999996E-03
0.500000 0.300000 4.999997E-03
0.600000 0.300000 4.999997E-03
0.700000 0.300000 4.999997E-03
0.800000 0.300000 4.999997E-03
0.900000 0.300000 4.999997E-03
-0.900000 0.400000 3.333333E-03
-0.800000 0.400000 4.444445E-03
-0.700000 0.400000 4.814815E-03
-0.600000 0.400000 4.938271E-03
-0.500000 0.400000 4.979424E-03
-0.400000 0.400000 4.993141E-03
-0.300000 0.400000 4.997713E-03
-0.200000 0.400000 4.999237E-03
-9.999999E-02 0.400000 4.999745E-03
1.490116E-08 0.400000 4.999915E-03
0.100000 0.400000 4.999971E-03
0.200000 0.400000 4.999990E-03
0.300000 0.400000 4.999996E-03
0.400000 0.400000 4.999998E-03
0.500000 0.400000 4.999999E-03
0.600000 0.400000 4.999999E-03
0.700000 0.400000 4.999999E-03
0.800000 0.400000 4.999999E-03
0.900000 0.400000 4.999999E-03
-0.900000 0.500000 3.333333E-03
-0.800000 0.500000 4.444445E-03
-0.700000 0.500000 4.814815E-03
-0.600000 0.500000 4.938272E-03
-0.500000 0.500000 4.979424E-03
-0.400000 0.500000 4.993142E-03
-0.300000 0.500000 4.997714E-03
-0.200000 0.500000 4.999238E-03
-9.999999E-02 0.500000 4.999746E-03
1.490116E-08 0.500000 4.999915E-03
0.100000 0.500000 4.999971E-03
0.200000 0.500000 4.999991E-03
0.300000 0.500000 4.999997E-03
0.400000 0.500000 4.999999E-03
0.500000 0.500000 4.999999E-03
0.600000 0.500000 5.000000E-03
0.700000 0.500000 5.000000E-03
0.800000 0.500000 5.000000E-03
0.900000 0.500000 5.000000E-03
-0.900000 0.600000 3.333333E-03
-0.800000 0.600000 4.444445E-03
-0.700000 0.600000 4.814815E-03
-0.600000 0.600000 4.938272E-03
-0.500000 0.600000 4.979424E-03
-0.400000 0.600000 4.993142E-03
-0.300000 0.600000 4.997714E-03
-0.200000 0.600000 4.999238E-03
-9.999999E-02 0.600000 4.999746E-03
1.490116E-08 0.600000 4.999916E-03
0.100000 0.600000 4.999972E-03
0.200000 0.600000 4.999991E-03
0.300000 0.600000 4.999997E-03
0.400000 0.600000 4.999999E-03
0.500000 0.600000 5.000000E-03
0.600000 0.600000 5.000000E-03
0.700000 0.600000 5.000000E-03
0.800000 0.600000 5.000000E-03
0.900000 0.600000 5.000000E-03
-0.900000 0.700000 3.333333E-03
-0.800000 0.700000 4.444445E-03
-0.700000 0.700000 4.814815E-03
-0.600000 0.700000 4.938272E-03
-0.500000 0.700000 4.979424E-03
-0.400000 0.700000 4.993142E-03
-0.300000 0.700000 4.997714E-03
-0.200000 0.700000 4.999238E-03
-9.999999E-02 0.700000 4.999746E-03
1.490116E-08 0.700000 4.999916E-03
0.100000 0.700000 4.999972E-03
0.200000 0.700000 4.999991E-03
0.300000 0.700000 4.999997E-03
0.400000 0.700000 4.999999E-03
0.500000 0.700000 5.000000E-03
0.600000 0.700000 5.000000E-03
0.700000 0.700000 5.000000E-03
0.800000 0.700000 5.000000E-03
0.900000 0.700000 5.000000E-03
-0.900000 0.800000 3.333333E-03
-0.800000 0.800000 4.444445E-03
-0.700000 0.800000 4.814815E-03
-0.600000 0.800000 4.938272E-03
-0.500000 0.800000 4.979424E-03
-0.400000 0.800000 4.993142E-03
-0.300000 0.800000 4.997714E-03
-0.200000 0.800000 4.999238E-03
-9.999999E-02 0.800000 4.999746E-03
1.490116E-08 0.800000 4.999916E-03
0.100000 0.800000 4.999972E-03
0.200000 0.800000 4.999991E-03
0.300000 0.800000 4.999997E-03
0.400000 0.800000 4.999999E-03
0.500000 0.800000 5.000000E-03
0.600000 0.800000 5.000000E-03
0.700000 0.800000 5.000000E-03
0.800000 0.800000 5.000000E-03
0.900000 0.800000 5.000000E-03
-0.900000 0.900000 3.333333E-03
-0.800000 0.900000 4.444445E-03
-0.700000 0.900000 4.814815E-03
-0.600000 0.900000 4.938272E-03
-0.500000 0.900000 4.979424E-03
-0.400000 0.900000 4.993142E-03
-0.300000 0.900000 4.997714E-03
-0.200000 0.900000 4.999238E-03
-9.999999E-02 0.900000 4.999746E-03
1.490116E-08 0.900000 4.999916E-03
0.100000 0.900000 4.999972E-03
0.200000 0.900000 4.999991E-03
0.300000 0.900000 4.999997E-03
0.400000 0.900000 4.999999E-03
0.500000 0.900000 5.000000E-03
0.600000 0.900000 5.000000E-03
0.700000 0.900000 5.000000E-03
0.800000 0.900000 5.000000E-03
0.900000 0.900000 5.000000E-03
1.00000 1.00000 0.00000 ]
GOING BACK TO QS, HOW DO I Set up the matrix system Au=b (taking the natural ordering of the nodes) and program the Jacobi iteration method to solve it, i.e.
u^k = (I-D^-1 A)u^(k-1) + D^-1 b, where k is the iteration index and D is the diagonal of A?
HOW DO I KNOW WHERE TO STOP AND THE ALGORITHM CONVERGES?
RE: Jacobi iterative method to solve elliptic equation
CODE
Refereing back to question, i think i have answer the question now with the code, but with this i have got some runtime error.
Run time error: Undefeined variables on arrays or element function :u_old(i,j)= u(i,j)
François Jacq : Dear mentor, iam trying to answer your questions as much as i can and i hope you have understand my code better this time. If not, please accept my apology. I understand what iam doing generally, but i am getting hard time to use it in fortran programming.
Note : There were no complier errors...
RE: Jacobi iterative method to solve elliptic equation
Anyway, your last version is very close to the solution except that you don't use explicitly the diagonal of the matrix as asked by your teacher.
I continue to be disturbed by the way you dimension your arrays. Compare u(0:ji+1,0:ji+1) and x(0:ji) for instance. You should have for any j, a X value for each U value. This is not dangerous here because x is in fact never used but...
If you continue to dimension u(0:ji+1,0:ji+1), which is a valid possible choice, then you must notice that you have ji+2 values in each direction and therefore ji+1 intervals. Because of the boundary conditions which are imposed, this leads to ji unknowns in each direction.
for a square -1. <= x <= 1. and -1. <= y <= 1., then the step size is h=2./(ji+1)
The array X should be declared (0:ji+1) and initialized as follows :
CODE --> fortran
Having now two arrays u and v is correct ! But I don't understand the reason of the third array u_old. You might also define the function you want to nullify; i.e. 1+L(u) (in fact this is not necessary for the Jacobi algorithm). Let us call F the corresponding array. As it is only valid inside the domain, and not at the boundary, its adjusted size is F(1:ji,1:ji). It could be constructed as follows :
CODE --> fortran
F is an affine transformation of u and may be interpreted as A.u-b. The matrix A is rather complicated. Fortunately, you don't have to construct it. The diagonal has the same size as F. You might declare it D(1:ji,1:ji). It just corresponds to dFij/duij => Dij=-4/h**2. As all the terms are constant, you don't need to construct the array D. We will see also that we don't have to construct F because the Jacobi algorithm is simple.
The Jacobi algorithm is mathematically : v = u - D^-1 F (The Newton algorithm would be v=u-A^-1 F). If one tries to simplify the relation, then one sees effectively that v(i,j) matches what you programmed ! Are you lucky ?
But to stay close to the text of your problem, you could write :
CODE
This is exactly equivalent to your version which is simpler :
CODE
François Jacq
RE: Jacobi iterative method to solve elliptic equation
CODE
I have no compliers error but i got runtime errors saying out of bounds on this lines:
u(i+1,j)=0 !East u(i+1,j), x direction
v(i,j)=0.25*(u(i+1,j)+u(i,j+1)+u(i-1,j)+u(i,j-1)+h**2)
u(i,j)=v(i,j) ! if do u=v, its say noncomfortable arrays
RE: Jacobi iterative method to solve elliptic equation
CODE
About the initialization of u, instead of using the loops which boundaries are wrong (the loop over i should go from 1 to ji), you should write without any loop :
CODE
About v(i,j), the wrong instruction is obviously u(i,j)=v(i,j) (written outside the loops over i an j => i==ji+1 and j==ji+1). First of all, you could allocate v exactly like u to make possible u=v. But even without changing anything else, the following instruction is ok :
CODE
For the final loops to print out the result, you should use do i=1,ji (or even do i=0,ji+1) and do j=1,ji (or even do j=0,ji+1) : why ji-1 for the upper bound ???
François Jacq
RE: Jacobi iterative method to solve elliptic equation
CODE -->
There are no errors and i can generate the outputs but i want to plot the graph of u(i,j) between -1<x<1 and -1<y<1 in matlab. My graph is not plotting properly between that axis especially in y-axis and the value of u(i,j) is not coming up properly on the graph.
RE: Jacobi iterative method to solve elliptic equation
CODE
This is the outputs when i run the program, when h=0.1 And iam trying to plot this using surf in the matlab but something is not working.
RE: Jacobi iterative method to solve elliptic equation
However, it seems not difficult to understand that you have to loop either over all the nodes (including the boundary) or over the nodes inside the domain (excluding the boundary).
In addition, you should take the time to improve the presentation of your program, especially the indentation of 2 or 3 characters to exhibit the structure, the enclosed loops... It is a mess to read it !
François Jacq
RE: Jacobi iterative method to solve elliptic equation