User:Yongtao-Li/FE-HW2.9

=Problem 2.9: "Example of WRF"

Given
General one dimensional Model 1.0/Data set 2 on P12-1:

where $$ {a_2} = 2,f = 3 $$

Consider $$ {b_j}\left( x \right) = \cos \left( {jx + \phi } \right) $$ and $$ \phi = \frac{1}{4}\pi $$ or $$ \phi  = \frac{3}{4}\pi $$

Find
A) Let n=2, What is the ndof of this problem?

B) Find 2 equations that enforce boundary conditions for $$ {u^h}\left( x \right) = \sum\limits_{j = 0}^n $$

C) Find 1 more equation to solve for $$ {\mathbf{d}} = {\left\{ \right\}_{3 \times 1}} $$ by projecting the residue $$ {\mathbf{P}}\left(  \right) $$ on a basis function $$ {b_k}(x)   with  k = 0,1,2  $$ s.t. the additional equation is linear independence from the above 2 equations

D) Display 3 equations in matrix form $$ {\mathbf{Kd}} = {\mathbf{F}} $$ and observe symmetric property of $$ {\mathbf{K}} $$

E) Solve for $$ {\mathbf{d}} $$

F) Construct $$ u_n^h\left( x \right) $$ and plot $$ u_n^h\left( x \right) \ vs. \ u\left( x \right) $$ ($$ u\left( x \right) $$ is the exact solution)

G) Repeat (1) to (6) for n=4 and n=6

H) Compare $$ u_n^h\left( {x = 0.5} \right)forn = 2,4,6 $$. If error $$ {e_n}\left( {0.5} \right) = u\left( {0.5} \right) - {u^h}\left( {0.5} \right) $$, plot $$ {e_n}\left( {0.5} \right) \ vs. \ n $$

Solution
A) Because n begins from 0, so the degree of freedom is 3 when n=2.

B) From (9.2) we can get $$ \sum\limits_{j = 0}^2 {{d_j}{b_j}\left( 1 \right)} = 0 $$

Then rewrite the above equation in matrix form:

where

From (9.3) we can get $$ \sum\limits_{j = 0}^2 {{d_j}b_j^{\left( 1 \right)}\left( 0 \right)} =  - 4 $$

Then rewrite the above equation in matrix form:

where

C) Because the residue $$ P = \frac{d}\left( {{a_2}\frac{d}{b_j}\left( x \right)} \right) - f\left( x \right)$$, after projecting the residue we get:

$$ \sum\limits_{j = 0}^2 {\left\{ {\int_0^1 {{b_k}\left( x \right)\left[ {\frac{d}\left( {{a_2}\frac{d}{b_j}\left( x \right)} \right)} \right]dx} } \right\}} {d_j} = - \int_0^1 {{b_k}\left( x \right)f\left( x \right)dx} $$

Then rewrite the above equation in matrix form:

where

E)-H) In order to simplify the computation procedure, the following Fortran code based on (9.8), (9.9) and (9.10) was used to calculate the results, when n=2,4,6 and $$\phi = \frac{1}{4}\pi,\frac{3}{4}\pi$$:

program main use lib implicit none integer :: n,i,j,k,m,l real,allocatable :: KK,d,F,b real :: fi(2),pi,x(1001),y(1001),yh(1001),e(3,2) n=6; pi=3.1415926; fi(1)=pi/4.0; fi(2)=pi*3.0/4.0 open (unit=11, file="egm6611.s11.team6.hw2.9.1.plt", status="replace") open (unit=12, file="egm6611.s11.team6.hw2.9.2.plt", status="replace") open (unit=13, file="egm6611.s11.team6.hw2.9.3.plt", status="replace") do l=2,n,2 write(*,*) "When the degree of freedom is ", l+1 do i=1,2 allocate(KK(l+1,l+1),d(l+1),F(l+1),b(2,l+1)) do j=0,l b(1,j+1)=cos(j*1.0+fi(i)) b(2,j+1)=-1.0*j*sin(j*0.0+fi(i)) end do     do k=0,l do j=0,l if (k==j) then KK(k+1,j+1)=-1.0*k/2.0*(2.0*k-sin(2*fi(i))+sin(2*fi(i)+2*k)) else KK(k+1,j+1)=2*j*j*j*(cos(fi(i))*sin(fi(i))-cos(fi(i)+k)*sin(fi(i)+j))-k*(cos(fi(i))*sin(fi(i))-cos(fi(i)+j)*sin(fi(i)+k)))/(j*j-k*k)           end if            if (k==0) then               F(k+1)=-3.0*cos(fi(i))            else               F(k+1)=-3.0*(sin(fi(i)+k)-sin(fi(i)))/k            end if         end do      end do      write(*,*) "when phi= ",fi(i)      do j=1,l+1      KK(1,j)=b(1,j)      KK(2,j)=b(2,j)      end do      F(1)=0     F(2)=-4     call invert(KK)     do k=1,l+1        d(k)=0.0        do j=1,l+1           d(k)=d(k)+KK(K,j)*F(j)        end do     end do     write(*,*) "solution vector d"     write(*,*) d     x=0.0; y=0.0; yh=0.0     do m=1,1001     x(m)=(m-1)/1000.0     y(m)=-3.0/4.0*x(m)*x(m)-4.0*x(m)+19.0/4.0     do k=1,l+1        yh(m)=yh(m)+d(k)*cos((k-1)*x(m)+fi(i))     end do   end do   do k=1,l+1      e(l/2,i)=e(l/2,i)+d(k)*cos((k-1)*0.5+fi(i)) end do  e(l/2,i)=(-3.0/4.0*0.5*0.5-4.0*0.5+19.0/4.0)-e(l/2,i) if (i==1) then write(11,*) "zone" do m=1,1001 write(11,"(3F12.4)") x(m),y(m),yh(m) end do  else write(12,*) "zone" do m=1,1001 write(12,"(3F12.4)") x(m),y(m),yh(m) end do  end if   deallocate(KK,d,F,b) end do !for i=1-2 end do !for n=1-8 write(13,"(A,2F12.4)") "2",e(1,1),e(1,2) write(13,"(A,2F12.4)") "4",e(2,1),e(2,2) write(13,"(A,2F12.4)") "6",e(3,1),e(3,2) end program

module lib implicit none contains !  subroutine invert(matrix) ! invert a small square matrix onto itself implicit none real,intent(in out)::matrix integer::i,k,n; real::con ; n= ubound(matrix,1) do k=1,n con=matrix(k,k); matrix(k,k)=1. matrix(k,:)=matrix(k,:)/con do i=1,n if(i/=k) then con=matrix(i,k); matrix(i,k)=0.0 matrix(i,:)=matrix(i,:) - matrix(k,:)*con end if     end do   end do   return end subroutine invert end module lib

Run this code and then get the output result for all the situations as following:

When the degree of freedom is           3 when phi=  0.7853981 solution vector d     1.142170       7.209867     -0.7765067 when phi=   2.356194 solution vector d    -17.24394       13.96153      -4.152339 When the degree of freedom is           5 when phi=  0.7853981 solution vector d    -2.044705       13.62609      -7.081253       2.811028     -0.5599513 when phi=   2.356194 solution vector d    -16.73928       16.00889      -7.965427       2.459494     -0.4499143 When the degree of freedom is           7 when phi=  0.7853981 solution vector d    -2.067435       14.06402      -8.531017       4.714233      -1.883228 0.4795609    -5.8918476E-02 when phi=   2.356194 solution vector d    -19.44370       26.57353      -25.56540       18.85248      -9.806686 3.187393    -0.5081997

For n=2,4,6, $$ \phi = \frac{1}{4}\pi $$, we compare the results with exact solution in the following image:



For n=2,4,6, $$ \phi = \frac{3}{4}\pi $$, we compare the results with exact solution in the following image:



In both the images, we can observe that the results become more close compared to exact solution when n increases.

For different shift phases, we compare the how the error's convergence situation in the following image:



It seems that both the shift phase can make a convergence around n=6.