第一章：线性方程组的解法 一、矩阵分解与线性方程组直接方法1. LU分解解方程module LU!module coment! Version : V1.0! Coded by : syz! Date :! Description : LU分解解方程!containssubroutine solve(A,b,x,N)implicit real*8(a-z)integer:Nreal*8:A(N,N),b(N),x(N)real*8:L(N,N),U(N,N)real*8:y(N)call doolittle(A,L,U,N)call downtri(L,b,y,N)call uptri(U,y,x,N)end subroutine solvesubroutine doolittle(A,L,U,N) !subroutine comment! Version : V1.0! Coded by : syz! Date :! Purpose : LU 分解之 Doolittle 函数!A=LU! Input parameters :!1.A方阵!2.N阶数! Output parameters :!1.L!2.U! Common parameters :! Post Script :!1.!2.!implicit real*8(a-z)integer:N,i,k,rreal*8:A(N,N),L(N,N),U(N,N)!U的第一行U(1,:)=A(1,:)此的第一列L(:,1)=a(:,1)/U(1,1)do k=2,Nl(k,k)=1do j=k,ns=0do m=1,k-1s=s+l(k,m)*u(m,j)end dou(k,j)=a(k,j)-send dodo i=k+1,ns=0do m=1,k-1s=s+l(i,m)*u(m,k)end dol(i,k)=(a(i,k)-s)/u(k,k)end doend subroutine doolittlesubroutine uptri(A,b,x,N)!subroutine comment! Version : V1.0! Coded by : syz! Date : 2010-4-8! Purpose :上三角方程组的回带方法!Ax=b! Input parameters :!1.A(N,N)系数矩阵!2.b(N)右向量!3.N方程维数! Output parameters :!1. x方程的根!2.! Common parameters :!implicit real*8(a-z)integer:i,j,k,Nreal*8:A(N,N),b(N),x(N)x(N)=b(N)/A(N,N)!回带部分do i=n-1,1,-1x(i)=b(i)do j=i+1,Nx(i)=x(i)-a(i,j)*x(j)end dox(i)=x(i)/A(i,i)end do end subroutine uptrisubroutine downtri(A,b,x,N)!subroutine comment! Version : V1.0! Coded by : syz! Date : 2010-4-9! Purpose :下三角方程组的回带方法!Ax=b! Input parameters :!1.A(N,N)系数矩阵!2.b(N)右向量!3.N方程维数! Output parameters :!1. x方程的根!2.! Common parameters :!implicit real*8(a-z)integer:i,j,Nreal*8:A(N,N),b(N),x(N)x(1)=b(1)/a(1,1)do k=2,Nx(k)=b(k)do i=1,k-1x(k)=x(k)-a(k,i)*x(i)end dox(k)=x(k)/a(k,k)end doend subroutine downtriend module LU program main!program comment! Version : V1.0! Coded by : syz! Date : 2010-4-8! Purpose : LU分解计算线性方程组!Ax=b! In put data files :!1.A,b!2.! Output data files :!1.x!2.!use LUinteger,parameter:N=4real*8:A(n,n),L(N,N),U(N,N)real*8:b(N),x(N)open(unit=11,file=fin.txt)open(unit=12,file=fout.txt)read(11,*)read(11,*)(A(i,j),j=1,N),i=1,N)read(11,*)bcall solve(A,b,x,N)write(12,101)x101 format(T5,LU分解计算线性方程组计算结果,/,4(/,F10.6)end program main2. LU分解之Croutmodule croutcontainssubroutine solve(A,L,U,N)!subroutine comment! Version : V1.0! Coded by : syz! Date :! Purpose : LU 之 Crout 分解!A=LU! Inputparameters :!1.A方阵!2.N阶数! Output parameters :!1.L!2.U! Common parameters :! Post Script :!1.!2.!implicit real*8(a-z)integer:N,i,k,rreal*8:A(N,N),L(N,N),U(N,N)!L第一列L(:,1)=a(:,1)!U第一行U(1,:)=a(1,:)/L(1,1) do k=2,Ndo i=k,ns=0do r=1,k-1s=s+l(i,r)*u(r,k)end dol(i,k)=a(i,k)-send dodo j=k+1,ns=0do r=1,k-1s=s+l(k,r)*u(r,j)end dou(k,j)=(a(k,j)-s)/l(k,k)end do u(k,k)=1end do end subroutine solveend module crout program main!program comment! Version : V1.0! Coded by : syz! Date : 2010-4-8! Purpose : Crot 分解! In put data files :!1.A,N!2.! Output data files :!1.L,U!2.!use croutinteger,parameter:N=4real*8:A(n,n),L(N,N),U(N,N)open(unit=11,file=fin.txt)open(unit=12,file=fout.txt)read(11,*)read(11,*)(A(i,j),j=1,N),i=1,N)call solve(A,L,U,N)write(12,21)21 format(T10,LU 之 Crout 分解,/)!输出L矩阵write(12,*)L=do i=1,Nwrite(12,22)L(i,:)end do22 format(4F10.6)!输出U矩阵write(12,*)U=do i=1,Nwrite(12,22)U(i,:)end do23 format(4F10.6)end program main3. LU 分解之 Doolittlemodule doolittle!module coment! Version : V1.0! Coded by : syz! Date :! Description : LU 分解之 doolittle 模块!containssubroutine solve(A,L,U,N)!subroutine comment! Version : V1.0! Coded by : syz! Date :! Purpose : LU 分解之 Doolittle 函数!A=LU! Inputparameters :!1.A方阵!2.N阶数! Output parameters :!1.L!2.U! Common parameters :! Post Script :!1.!2.!implicit real*8(a-z)integer:N,i,k,rreal*8:A(N,N),L(N,N),U(N,N)!U的第一行U(1,:)=A(1,:)此的第一列L(:,1)=a(:,1)/U(1,1)do k=2,Nl(k,k)=1do j=k,ns=0do m=1,k-1s=s+l(k,m)*u(m,j)end dou(k,j)=a(k,j)-send dodo i=k+1,ns=0do m=1,k-1s=s+l(i,m)*u(m,k)end dol(i,k)=(a(i,k)-s)/u(k,k)end do end doend subroutine solveend module doolittle program main!program comment! Version : V1.0! Coded by : syz! Date : 2010-4-8! Purpose : Doolittle 分解! In put data files :!1.A,N!2.! Output data files :!1.L,U!2.!use doolittleinteger,parameter:N=3real*8:A(n,n),L(N,N),U(N,N)open(unit=11,file=fin.txt)open(unit=12,file=fout.txt)read(11,*)read(11,*)(A(i,j),j=1,N),i=1,N)call solve