代码拉取完成,页面将自动刷新
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);
if(r>0.5)
disp('r>0.5,unstability');
end
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
% for n = 1:N + 1
% u([1 M],n) = [bx0(t(n)); bxf(t(n))];
% end
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)
% for l=2:M
% if(u(l,k)>u(l-1,k))
% u(l,k)=u(l-1,k);
% end
% end
end
u=u';
此处可能存在不合适展示的内容,页面不予展示。您可通过相关编辑功能自查并修改。
如您确认内容无涉及 不当用语 / 纯广告导流 / 暴力 / 低俗色情 / 侵权 / 盗版 / 虚假 / 无价值内容或违法国家有关法律法规的内容,可点击提交进行申诉,我们将尽快为您处理。