二阶椭圆偏微分方程实例求解(附matlab代码)(2)

for i=1:2*n-1 %使A.B下标i-1/2变为2i-1 for j=1:2*n-1 x(i)=i*h/2; y(j)=j*h/2;

A(i,j)=FA(x(i),y(j)); B(i,j)=FB(x(i),y(j)); end end

x=zeros(n-1,1); y=zeros(n-1,1); for i=1:n-1 for j=1:n-1 x(i)=i*h; y(j)=j*h;

C(i,j)=FC(x(i),y(j)); D(i,j)=FD(x(i),y(j)); E(i,j)=1;

F(i,j)=FF(x(i),y(j)); U(i,j)=FU(x(i),y(j)); end end

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

%对最终要解的方程组的系数矩阵进行赋值 for i=1:n-1

bcd=zeros(n-1); bb=[]; cc=[]; dd=[]; gg=[]; kk=[]; for j=1:n-1

b(i,j)=(A(2*i+1,2*j)+A(2*i-1,2*j)+B(2*i,2*j+1)+B(2*i,2*j-1))/h^2+E(i,j); c(i,j)=(B(2*i,2*j+1)-h*D(i,j)/2)/h^2;

d(i,j)=(B(2*i,2*j-1)+h*D(i,j)/2)/h^2; g(i,j)=(A(2*i+1,2*j)-h*C(i,j)/2)/h^2;

k(i,j)=(A(2*i-1,2*j)+h*C(i,j)/2)/h^2; bb=[bb b(i,j)]; if j<=n-2

cc=[cc c(i,j)]; end if j>=2

dd=[dd d(i,j)]; end

gg=[gg g(i,j)]; kk=[kk k(i,j)];

%给f赋值 if i==1

ff(i,j)=F(i,j)+k(i,j)*1;%边值为1 elseif i==n-1

ff(i,j)=F(i,j)+g(i,j)*A(i,2*j);%A中i取值无所谓，不影响 else

ff(i,j)=F(i,j); end if j==1

ff(i,j)=ff(i,j)+d(i,j)*1;%边值为1 elseif j==n-1

ff(i,j)=ff(i,j)+c(i,j)*B(2*i,j); end

f((i-1)*(n-1)+j,1) = ff(i,j);%你应该懂的，坐标变换 end

bcd=diag(bb)-diag(cc,1)-diag(dd,-1); BCD=blkdiag(BCD,bcd); if i<=n-2

G=[G gg]; end

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

if i>=2

K=[K kk]; end end

P=BCD-diag(G,n-1)-diag(K,-n+1);

%BCD=BCD-diag(G,n-1)-diag(K,-n+1);x=Doolittle(BCD,f);这样是不是可以减少点内存消耗

x=Doolittle(P,f);%LU分解解方程组，这里也可以输入其他解方程组的方法 %x=GaussXQByOrder(P,f) %高斯消去法

%x0=ones((n-1)^2,1);x=gauseidel(P,f,x0) %gauseidel迭代法 for i=1:n-1 for j=1:n-1

u((i-1)*(n-1)+j)=U(i,j); end end

error=abs(x-u); result=[x u error]; time = toc;

disp('计算时间为：'); disp(time);

disp('---------------------------------------------------------------'); format long;

disp('计算结果为：');

disp(' 数值解真实值误差'); disp(result);

6.2LU分解解线性方程组

function [x,L,U]= Doolittle (A,b) N = size(A); n = N(1);

L = eye(n,n); %L的对角元素为1 U = zeros(n,n);

U(1,1:n) = A(1,1:n); %U的第一行 L(1:n,1) = A(1:n,1)/U(1,1); %L的第一列

for k=2:n for i=k:n

U(k,i) = A(k,i)-L(k,1:(k-1))*U(1:(k-1),i); %U的第k行 end

for j=(k+1):n

L(j,k) = (A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k);

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

%L的第k列 end end

y = SolveDownTriangle(L,b);

x = SolveUpTriangle(U,y); %求解方程

6.3高斯消去法

function [x,XA]=GaussXQByOrder(A,b) N = size(A); n = N(1);

for i=1:(n-1) for j=(i+1):n if(A(i,i)==0)

disp('对角元素为0！'); %防止对角元素为0 return; end

l = A(j,i); m = A(i,i);

A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m; %消元方程 b(j)=b(j)-l*b(i)/m; end end

x=SolveUpTriangle(A,b); %通用的求上三角系数矩阵线性方程组的函数 XA = A; %消元后的系数矩阵

6.4上三角解线性方程组（LU分解法、高斯消去法要用到这个算法） function x=SolveUpTriangle(A,b) N=size(A); n=N(1); for i=n:-1:1 if(i<n)

s=A(i,(i+1):n)*x((i+1):n,1); else s=0; end

x(i,1)=(b(i)-s)/A(i,i); end

6.5下三角解线性方程组（LU分解法要用到这个算法） function x=SolveDownTriangle(A,b)

二阶椭圆偏微分方程实例求解(附matlab代码)用的是五点差分法。

N=size(A); n=N(1); for i=1:n if(i>1)

s=A(i,1:(i-1))*x(1:(i-1),1); else s=0; end

x(i,1)=(b(i)-s)/A(i,i); end

6.6gauseidel迭代法 ifnargin==3

eps= 1.0e-6;%可以安自己的精度要求改变 M =10000;

elseifnargin == 4 M =10000;

elseifnargin<3 error return; end

D=diag(diag(A)); %求A的对角矩阵 L=-tril(A,-1); %求A的下三角阵 U=-triu(A,1); %求A的上三角阵 G=(D-L)\U; f=(D-L)\b; x=G*x0+f;

n=1; %迭代次数 while norm(x-x0)>=eps x0=x;

x=G*x0+f; n=n+1; if(n>=M)

disp('Warning: 迭代次数太多，可能不收敛！'); return; end end

disp('迭代次数为：'); disp(n);