Tham khảo tài liệu 'shallow liquid simulation using matlab (2001 neumann) episode 4', kỹ thuật - công nghệ, cơ khí - chế tạo máy phục vụ nhu cầu học tập, nghiên cứu và làm việc hiệu quả | ------------------------------------------------------- Given vorticity omega Solve for streamfunction psi. if solver use_fft if extra_row matA 1 1 matA 1 1 cludge end end if solver use_fft Take FFT of omega in x then in y. note fft works along columns so transpose afterwards to work on other dimension and transpose back at end. ft fft fft omega . . Then divide by kx 2 ky 2 where kx ky are x y values on the grid in frequency space . ft ft. k1 Then take iFFT in x then in y to get psi. psi ifft ifft ft . . psi real psi elseif solver use_backslash solve using backslash psi matA reshape omega N2 1 elseif solver use_lu if extra_row y1 matL reshape omega N2 1 0 else y1 matL reshape omega N2 1 end psi matU y1 elseif solver use_bicg pass current value of psi as estimate for next iteration psi bicgstab matA reshape omega N2 1 1e-4 1000 . reshape psi N2 1 elseif solver use_gmres pass current value of psi as estimate for next iteration psi gmres matA reshape omega N2 1 1e-4 1000 . reshape psi N2 1 end if solver use_fft if extra_row matA 1 1 matA 1 1 - cludge end psi reshape psi n n end 31 ------------------------------------------------------- Check whether del 2 psi omega Use central difference formula for 2nd derivative. if psi_2nd_diff psi_xx psi_yy omega So compute numerical 2nd deriv in x y directions and add. Result should match omega. save_n n n size psi 1 use central difference for 2nd derivatives df psi 1 n-2 2 n-1 -2. psi 2 n-1 2 n-1 psi 3 n 2 n-1 df df psi 2 n-1 1 n-2 -2. psi 2 n-1 2 n-1 psi 2 n-1 3 n df df dX dx n save_n df df - omega 2 n-1 2 n-1 fprintf diff at omega 1 1 e n df 1 1 df df - df 1 1 figure 1 pcolor seems to not show the last row column df2 df df 1 df2 df2 df2 1 pcolor df2 colorbar shading flat drawnow max max df min min df clear save_n df psi2 fprintf break n break end ------------------------------------------------------ set up matrix B corresponds to the curl product of psi cross omega strategy 1. create a B matrix from loops. big slow but .
