急求！！！ 大学数学，用matlab解决问题，题目是一维抛物型偏微分方程差分解法
1,不一定有效果，因为pdetool具体编程是不知道的，如果解决小问题两者的结果一样说明不了什麽问题，尤其对于偏微分方程。
2有限元的边界必须固定，从数理方程上讲静态有限元问题就是边值问题，如果边界变化的话，初始一下别的专业有限元软件，比如anasys，adima等。
MATLAB提供两种解决PDE问题:
pdepe()函数求解般PDEs据用较通用性支持命令行形式调用
二PDE工具箱求解特殊PDE问题PDEtool较局限性比能求解二阶PDE问题并且能解决偏微程组提供GUI界面繁杂编程解脱同通File->Save As直接M代码
MATLAB语言提供pdepe()函数直接求解般偏微程(组)调用格式
sol=pdepe(m,@pdefun,@pdeic,@pdebc,x,t)
【输入参数】
@pdefun:PDE问题描述函数必须换面标准形式
PDE编写面入口函数
[c,f,s]=pdefun(x,t,u,du)
m,x,t应于(式1)相关参数duu阶导数由给定输入变量即表示c,f,s三函数
@pdebc:PDE边界条件描述函数必须先化面形式
于边值条件编写面函数描述
[pa,qa,pb,qb]=pdebc(x,t,u,du)
其a表示边界b表示边界
@pdeic:PDE初值条件必须化面形式
股我使用面简单函数描述
u0=pdeic(x)
m,x,t:应于(式1)相关参数
【输参数】
sol:三维数组sol(:,:,i)表示ui解换句说uk应x(i)t(j)解sol(i,j,k)
通sol我使用pdeval()直接计算某点函数值
显式前向欧拉法
源程序：
function [u,x,t] = EF_Euler(A,xf,T,it0,bx0,bxf,M,N)
%解方程 A u_xx = u_t ， 0 <= x <= xf, 0 <= t <= T
%初值: u(x,0) = it0(x)
% 边界条件: u(0,t) = bx0(t), u(xf,t) = bxf(t)
% M ： x 轴的等分段数
% N ： t 轴的等分段数
dx = xf/M; x = [0:M]*dx;
dt = T/N; t = [0:N]'*dt;
for i= 1:M + 1
u(i,1) =it0(x(i));
end
for j = 1:N + 1
u([1 M + 1],j)= [bx0(t(j)); bxf(t(j))];
end
r = A*dt/dx/dx, r1 = 1 - 2*r;
if(r>0.5)
disp('r>0.5,unstability');
end
for j = 1:N
for i = 2:M
u(i,j+1) =r*(u(i + 1,j) + u(i-1,j)) + r1*u(i,j); %(9.2.3)
end
end
u=u';
%------------------------------------
在MATLAB中编写脚本文件：
A = 0.5; %方程系数
it0 = inline('sin(pi*x)','x'); %初始条件
bx0 = inline('0'); bxf = inline('0'); %边界条件
xf = 2; M = 80; T = 0.1; N = 100;
[u1,x,t] = EF_Euler(A,xf,T,it0,bx0,bxf,M,N);
figure(1),clf,mesh(u1)
xlabel('x')
ylabel('t')
zlabel('U')
title('r>0.5')
M=50;
[u1,x,t] = EF_Euler(A,xf,T,it0,bx0,bxf,M,N);
figure(2),clf,mesh(u1)
xlabel('x')
ylabel('t')
zlabel('U')
title('r<0.5')
%-----------------
隐式后向欧拉法
源程序：
function [u,x,t] = IB_Euler(A,xf,T,it0,bx0,bxf,M,N)
%解方程 A1 u_xx = u_t ， 0 <= x <= xf, 0 <= t <= T
%初值: u(x,0) = it0(x)
% 边界条件: u(0,t) = bx0(t), u(xf,t) = bxf(t)
% M ： x 轴的等分段数
% N ： t 轴的等分段数
dx = xf/M; x = [0:M]*dx;
dt = T/N; t = [0:N]'*dt;
for i = 1:M + 1
u(i,1) =it0(x(i));
end
for j = 1:N + 1
u([1 M + 1],j)= [bx0(t(j)); bxf(t(j))];
end
r = A*dt/dx/dx;
r2 = 1 + 2*r;
for i = 1:M - 1
P(i,i) = r2; %构造9.2.9的矩阵
if i > 1
P(i - 1,i)= -r; P(i,i - 1) = -r;
end
end
for j = 2:N + 1
b = [r*u(1,j);zeros(M - 3,1); r*u(M + 1,j)] + u(2:M,j - 1); %Eq.(9.2.9)
u(2:M,j) =linsolve(P,b);
end
u=u';
%---------------------
在MATLAB中编写脚本文件：
A= 0.5; %方程系数
it0 = inline('sin(pi*x)','x'); %初始条件
bx0 = inline('0'); bxf = inline('0'); %边界条件
xf = 2; M = 50; T = 0.1; N = 100;
[u1,x,t] = IB_Euler(A,xf,T,it0,bx0,bxf,M,N);
mesh(u1)
xlabel('x')
ylabel('t')
zlabel('U')
%--------------
Grank-Nicholson方法
源程序：
function [u,x,t] =Grank_Nicholson(A,xf,T,it0,bx0,bxf,M,N)
%解方程 A u_xx = u_t ， 0 <= x <= xf, 0 <= t <= T
%初值: u(x,0) = it0(x)
% 边界条件: u(0,t) = bx0(t), u(xf,t) = bxf(t)
% M ： x 轴的等分段数
% N ： t 轴的等分段数
dx = xf/M; x = [0:M]*dx;
dt = T/N; t = [0:N]'*dt;
for i = 1:M + 1
u(i,1) =it0(x(i));
end
for n = 1:N + 1
u([1 M + 1],n)= [bx0(t(n)); bxf(t(n))];
end
r = A*dt/dx/dx;
r1 = 2*(1 + r); r2 = 2*(1 - r);
for i = 1:M - 1
P(i,i) = r1; %(9.2.17)
Q(i,i)=r2;
if i > 1
P(i - 1,i)= -r; P(i,i - 1) = -r; %(9.2.17)等式左边矩阵
Q(i - 1,i)= r; Q(i,i - 1) = r;%(9.2.17)等式右边矩阵
end
end
for k = 2:N + 1
b=Q*u(2:M,k-1)+[r*(u(1,k)+u(1,k-1));zeros(M-2,1)];
u(2:M,k) =linsolve(P,b); %(9.2.17)
end
u=u';
%-----------
在MATLAB中编写脚本文件：
A= 0.5; %方程系数
it0 = inline('sin(pi*x)','x'); %初始条件
bx0 = inline('0'); bxf = inline('0'); %边界条件
xf = 2; M = 25; T = 0.1; N = 100;
[u1,x,t] = Grank_Nicholson(A,xf,T,it0,bx0,bxf,M,N);
mesh(u1)
xlabel('x')
ylabel('t')
zlabel('U')
急求！！！ 大学数学，用matlab解决问题，题目是一维抛物型偏微分方程差分解法视频
相关评论：
