matlab_偏微分方程数值解法中用C-N格式求解抛物型方程程序

%本程序用于Crank-Nicolson求解抛物型方程 %其中函数u|t=0 = g(x), u|x=0 = g1(t), u|x=1 = g2(t)在m文件中定义 function w=cn(M,h,k) %M = 10; h = 0.1; % M*h = 1 %k = 0.1; r = k/h^2; %去t方向步长0.1，同时t<=1 U = zeros(M,M-1); %用于存放所求的点 %向量u用于存放第0时间层的各点值 u = zeros(M-1,1); for i = 1:M-1 u(i) = g(i*h); end %给矩阵Ａ赋值 A = zeros(M-1,M-1); for i = 1:M-1 if i == 1 A(1,1) = 1 + r; A(1,2) = -0.5*r; else if i == (M-1) A(M-1,M-2) = -0.5*r; A(M-1,M-1) = 1 + r; else A(i,i-1) = -0.5*r; A(i,i) = 1 + r; A(i,i+1) = -0.5*r; end end end %给矩阵Ｂ赋值 B = zeros(M-1,M-1); for i = 1:M-1 if i == 1 B(1,1) = 1 - r; A(1,2) = 0.5*r; else if i ==(M-1) A(M-1,M-2) = 0.5*r; A(M-1,M-1) = 1 - r; else A(i,i-1) = 0.5*r; A(i,i) = 1 - r; A(i,i+1) = 0.5*r; end end end %求解个时间层的结点 for i = 1:M %求解第一时间层 if i==1 b = zeros(M-1,1); b(1) = 0.5*r*g1(0) + 0.5*r*g1(k); b(9) = 0.5*r*g2(0) + 0.5*r*g2(k); b = B*u + b; %用追赶法求解　A*x = b f = zeros(1,M-1); w = zeros(1,M-1); f(1) = b(1)/(1+r); w(1) = 0.5*r/(1+r); for j = 2:M-2 f(j) = (b(j) + 0.5*r*f(j-1))/(1+r - 0.5*r*w(j-1)); w(j) = 0.5*r/(1+r - 0.5*r*w(j-1)); end f(M-1) = (b(M-1) + 0.5*r*f(M-2))/(1+r - 0.5*r*w(M-2)); U(1,M-1) = f(M-1); for m = M-2:-1:1 U(1,m) = f(m) + w(m)*U(1,m+1); end else %求解其他时间层结点值 b = zeros(M-1,1); b(1) = 0.5*r*g1((i-1)*k) + 0.5*r*g1(i*k); b(M-1) = 0.5*r*g2((i-1)*k) + 0.5*r*g2(i*k); %vec为第i－1时间层的节点值 vec = zeros(M-1,1); for m = 1:M-1 vec(m) = U(i-1,m); end b = B*vec + b; %用追赶法求解　A*x = b f = zeros(1,M-1); w = zeros(1,M-2); f(1) = b(1)/(1+r); w(1) = 0.5*r/(1+r); for j = 2:M-2 f(j) = (b(j) + 0.5*r*f(j-1))/(1+r - 0.5*r*w(j-1)); w(j) = 0.5*r/(1+r - 0.5*r*w(j-1)); end f(M-1) = (b(M-1) + 0.5*r*f(M-2))/(1+r - 0.5*r*w(M-2)); U(i,M-1) = f(M-1); for m = M-2:-1:1 U(i,m) = f(m) + w(m)*U(i,m+1); end end end disp('输出各结点值'); U