ex Sverdrup equation Keywords Gauss elimination LU decomposition Steady equation We consider the steady linear equation If we dont need to know the temporal evolution it is better to directly solve the steady equation ID: 1046390
Download Presentation The PPT/PDF document "Solving linear equations" is the property of its rightful owner. Permission is granted to download and print the materials on this web site for personal, non-commercial use only, and to display it on your personal computer provided you do not modify the materials and that you retain all copyright notices contained in the materials. By downloading content from our website, you accept the terms of this agreement.
1. Solving linear equationsex) Sverdrup equationKeywordsGauss eliminationLU decomposition
2. Steady equationWe consider the steady linear equationIf we don’t need to know the temporal evolution, it is better to directly solve the steady equation By using the stream function , we obtainwhere2
3. Finite differenceWe consider the domain thatzonal: i = 1, 2, …, N (i=1 and N are the land)meridional: j = 1, 2, …, M (j=1 and M are the land)grid space: dx = dy The boundary condition is = 0 along the boundariesBy using these notations,3
4. Matrix formThese terms are changed into the matrix form4
5. Matrix formThen, the Sverdrup equation is represented bywhere A is a square matrix (M N), and and F are a vector [(M N)1]5unknown
6. Matrix formAlong the boundary, = 0The number of equations are 2N+2M‒4Except for the boundary, The number of equations are (N‒2)(M‒2)There are NM equations and NM unknown valuesWe can obtain the solution 6
7. Setting7Parameter settinggrid size: 100 kmnumber of x-grid: 48number of y-grid: 38 plain (f0 =10‒4 second‒1 and =10‒11 second‒1 m‒1)depth: 5000 m: 10‒6 second‒1Wind forcing x = ‒0.01cos(y/38) N/m2 y = 0
8. Example 1/3 in Matlabdx=100*1000; % mdy=dx;f0=1.0e-04; beta=1.0e-11;rho=1000;grev=9.8;alpha=1e-6;H=5000; % depth meterxsz=48; ysz=38;% initial conditionyy=0:dy:(dy*(ysz)); yy=yy./(dy*(ysz));taux=-0.01*ones(xsz,1)*cos(pi*yy);curltau=-(taux(:,2:end)-taux(:,1:end-1))./dy./rho./H;A=zeros(xsz*ysz,xsz,ysz); F=zeros(xsz*ysz,1);8
9. Example 2/3 in Matlabfor jj=1:ysz for ii=1:xsz aa=ii+(jj-1)*xsz; if(ii==1 || ii==xsz || jj==1 || jj==ysz) A(aa,ii,jj)=1; % boundary condition else A(aa,ii,jj)=-4.*alpha./dx./dx; A(aa,ii-1,jj)=alpha./dx./dx-beta./2./dx; A(aa,ii+1,jj)=alpha./dx./dx+beta./2./dx; A(aa,ii,jj-1)=alpha./dx./dx; A(aa,ii,jj+1)=alpha./dx./dx; F(aa,1)=curltau(ii,jj); end endend9
10. Example 3/3 in Matlab% solveA=reshape(A,xsz*ysz,xsz*ysz);phi=mldivide(A,F); % solve by Matlabphi=reshape(phi,xsz,ysz);h=f0.*phi./grev; % calculate thickness anomalyh=h-mean(h(:));10
11. Example 1/3 in Pythonimport numpy as npdx, dy=100.0*1000.0, 100.0*1000.0 # mf0, beta=1.0e-04, 1.0e-11rho=1000.0grev=9.8alpha=1.0e-6H=5000.0 # depth meterxsz, ysz=48, 38# initial conditiontaux = np.zeros((xsz, ysz+1)) for ii in range(0,xsz): for jj in range(0,ysz+1): taux[ii, jj] = -0.01*np.cos(np.pi*jj/(ysz-1))curltau=-(taux[:,1:]-taux[:,0:ysz])/dy/rho/H11
12. Example 2/3 in PythonA = np.zeros((xsz*ysz,xsz,ysz))F = np.zeros((xsz*ysz,1))for ii in range(0,xsz): for jj in range(0,ysz): aa = ii+jj*xsz if(ii==0 or ii==xsz-1 or jj==0 or jj==ysz-1): A[aa,ii,jj] = 1.0 # boundary condition else: A[aa,ii,jj] = -4.0*alpha/dx/dx A[aa,ii-1,jj] = alpha/dx/dx - beta/2.0/dx A[aa,ii+1,jj] = alpha/dx/dx + beta/2.0/dx A[aa,ii,jj-1] = alpha/dx/dx A[aa,ii,jj+1] = alpha/dx/dx F[aa,0] = curltau[ii,jj]12
13. Example 3/3 in PythonA=np.reshape(A,(xsz*ysz,xsz*ysz),order=“F”)invA=np.linalg.inv(A) # solve in Pythonphi=np.dot(invA,F)h=f0*phi/grev # calculate thickness anomalyh=h-np.mean(h)h=np.reshape(h,(xsz,ysz),order=“F”)13Here, we calculate
14. Result14xy[cm]thickness anomaly
15. How to solveThere are many methods to solve a set of linear equationsGauss elimination methodLU decompositionJacobi methodGauss-Seidel methodetc.We employ the Gauss elimination method, which is a basic method, and the LU decomposition, which is used in the Matlab “mldivide” function15
16. Gauss eliminationAs an example, we solve the equations belowTo solve these equations, they are transformed into the upper triangular matrix16
17. Gauss eliminationWe calculate (2) plus (1) and (3) minus 2(1) Next, we calculate (6) plus (5)17(1)(2)(3)(4)(5)(6)
18. Gauss eliminationFrom (9), z = 3From (8) and z = 3, y = 2From (7), y = 2 and z = 3, x = 118(7)(8)(9)Gauss elimination method
19. PivotHowever, if a diagonal component is zero...Before manipulating the equations, (10) is replaced to (12), before (12) has the maximum value in the first columnThis manipulation is called the pivotIn each step, it is better to do this procedure19
20. LU decompositionWe separate the matrix A into the lower triangular matrix L and the upper triangular matrix UIn this case, we can easily solve these two sets of the equationsA nice point is that L and U are unique to AThus, even if F is changed, we don’t need to calculate L and U again 20and
21. LU decompositionAs an example, we consider a 3 by 3 matrix, then L and U are as follows:Therefore,21
22. LU decompositionIn a more general case, at the i row,Since the denominator is lii, the pivot is sometimes necessary22j ≤ ij > i
23. Example 1/223First, we solve
24. Example 2/224Next, we solve
25. ClassificationPartial differential equationellipticPoisson’s equationLaplace’s equationSverdrup equationparabolicheat equationhyperbolicwave equation25Solve a set of linear equationsCalculate temporal integration of equations