ADI解二维抛物方程 matlab代码

最新推荐文章于 2024-04-15 21:24:50 发布

转载最新推荐文章于 2024-04-15 21:24:50 发布 · 3.6k 阅读

CC 4.0 BY-SA版权

原文链接：https://wenku.baidu.com/view/5c82173da32d7375a41780aa.html?fr=search

文章标签：

#线性代数 #算法 #矩阵

读代码专栏收录该内容

1 篇文章

订阅专栏

本文介绍使用交替方向隐式（ADI）方法求解二维热传导方程的具体步骤，包括数值解的迭代过程、误差分析及精确解的对比。通过MATLAB实现ADI算法，展示不同时间步长下的数值解与精确解，并进行图像绘制，最后分析两种解的误差。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

解以下方程
待解方程
用P-R ADI 格式
对于Laplace算子的处理

第一步对xx用显式，yy用隐式。
第二步对yy用显式，xx用隐式。

%**************************************************************
%MAIN.m
clc
clear
a=0; b=1;     %x取值范围
c=0; d=1;       %y取值范围
tfinal = 1;        %最终时刻
t=1/1600;%时间步长；
h=1/40;%空间步长
r=t/h^2;%网比
x=a:h:b;
y=c:h:d;
%精确解
m=40;
u1=zeros(m+1,m+1);
for i=1:m+1,
    for j=1:m+1
        u1(j,i) = uexact(x(i),y(j),1);
    end
end
%数值解
u=ADI(a,b,c,d,t,h,tfinal);
%绘制图像
figure(1) ;mesh(x,y,u1)
figure(2); mesh(x,y,u)
%误差分析
error=u-u1;
norm1=norm(error,1);
norm2=norm(error,2);
norm00=norm(error,inf);

%**************************************************************
%ADI.m
function [u]=ADI(a,b,c,d,t,h,tfinal )
 %(a , b) x取值范围
 %(c, d)  y取值范围
%tfinal最终时刻
%t时间步长；
%h空间步长
r=t/h^2;%网比
m=(b-a)/h;%
n=tfinal/t; %
x=a:h:b;
y=c:h:d;
%初始条件
u=zeros(m+1,m+1);
for i=1:m+1,
    for j=1:m+1
        u(j,i) = uexact(x(i),y(j),0);
    end
end
u2=zeros(m+1,m+1);
a=-1/32*r*ones(1,m-2);
b=(1+r/16)*ones(1,m-1); 
aa=-1/32*r*ones(1,m);
cc=aa;aa(m)=-1;cc(1)=-1;
bb=(1+r/16)*ones(1,m+1);
bb(1)=1;bb(m+1)=1;
for i=1:n % n time steps
    %从n->n+1/2,u_{xx}implicit scheme，u_{yy}explicit scheme
    for j=2:m
        for k=2:m
            d(k-1)=1/32*r*(u(j,k+1)-2*u(j,k)+u(j,k-1))+u(j,k);
        end
        % first and last term of d is zero.
        %d(1)=d(1)+u1(j,1);d(m-1)=d(m-1)+u1(j,m+1);
        u2(j,2:m)=Thomas(a,b,a,d);
    end
    u2(1,:)=u2(2,:);%u_y=0 at the boundary, if not assign value then 0
    u2(m+1,:)=u2(m,:);
    
    %从n->n+1,u_{xx}explicit，u_{yy}implicit
    for k=2:m %k=1, k=m+1, boundary condition=0, no need to update
        dd(1)=0;dd(m+1)=0;
        for j=2:m
            dd(j)=1/32*r*(u2(j+1,k)-2*u2(j,k)+u2(j-1,k))+u2(j,k);
        end
        u(:,k)=Thomas(aa,bb,cc,dd);
    end
    u2=u;
end
end

%**************************************************************
%Thomas.m
function [x]=Thomas(a,b,c,d)
%a: -1 diagonal, below diagonal line
%b: diagnal
%c: 1 diagonal, above diagonal line
%d: RHS, (diag(b)+diag(a,-1)+diag(b,1))x=d
r=size(a);
m=r(2);
r=size(b);
n=r(2);
if size(a)~=size(c)|m~=n-1|size(b)~=size(d)
    error('Incompatible input dimensions !');
end

%LU decomposition
%L=
%U=diag(ones())+diag(u,1)
u(1)=b(1);
for i=2:n
    l(i-1)=a(i-1)/u(i-1);
    u(i)=b(i)-l(i-1)*c(i-1);
    %v(i-1)=(b(i)-u(i))/l(i-1);   
end

%Ly=d
y(1)=d(1);
for i=2:n
    y(i)=d(i)-l(i-1)*y(i-1);
end

%Ux=y
x(n)=y(n)/u(n);
for i=n-1:-1:1
   x(i)=y(i)/u(i);
   x(i)=(y(i)-c(i)*x(i+1))/u(i);
end
end

%**************************************************************
%uexact.m
function [ f]=uexact(x,y,t)
f=sin(x*pi)*cos(y*pi)*exp(-pi*pi/8*t);
end