Fpg
һάÀ©É¢·½³Ì¤ÎÓÐÏÞ²î·Ö·¨
¡ª¡ª¼ÆËãÎïÀíʵÑé×÷ÒµÆß
³ÂÍò ÎïÀíѧ2013¼¶ 13020011006
? ÌâÄ¿£º
±à³ÌÇó½âһάÀ©É¢·½³Ì¤Î½â
??u?2u??D2(0?x?a0,0?t?tmax)?x??t?u(x,t)|t?0?ex???u au?b?c(x?0)?111?n??u?au?b?c2(x?a0)22??n?È¡a1?1,b1?1,c1?0,a2?1,b2??1,c2?0,a0?1.0,tmax?10,D?0.1,h?0.1,??0.1¡£Êä³öt=1,2,...,10ʱ¿Ì¤ÎxºÍu(x),²¢Óë½âÎö½âu=exp(x+0.1t)×÷±È½Ï¡£
? Ö÷³ÌÐò£º
% һάÀ©É¢·½³Ì¤ÎÓÐÏÞ²î·Ö·¨ clear,clc;
%¶¨Òå³õʼ³£Á¿
a1 = 1; b1 = 1; c1 = 0; a2 = 1;b2 = -1; c2 = 0; a0 = 1.0; t_max = 10; D = 0.1; h = 0.1; tao = 0.1;
%µ÷ÓÃÀ©É¢·½³Ì×Óº¯ÊýÇó½â
u = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2);
? ×Ó³ÌÐò1£º
function output = diffuse_equation(a0,t_max,h,tao,D,a1,b1,c1,a2,b2,c2)
% һάÀ©É¢·½³Ì¤ÎÓÐÏÞ²î·Ö·¨£¬²ÉÓÃÒþʽÁùµã²î·Ö¸ñʽ(Crank-Nicolson) % a0: x¤Î×î´óÖµ % t:_max: t¤Î×î´óÖµ % h: ¿Õ¼ä²½³¤ % tao: ʱ¼ä²½³¤ % D£ºÀ©É¢ÏµÊý
% a1,b1,c1ÊÇ£¨x=0£©±ß½çÌõ¼þ¤ÎϵÊý£»a2,b2,c2ÊÇ£¨x=a0£©±ß½çÌõ¼þ¤ÎϵÊý
x = 0:h:a0;
Fpg
Fpg
n = length(x); t = 0:tao:t_max; k = length(t);
P = tao * D/h^2; P1 = 1/P + 1; P2 = 1/P - 1;
u = zeros(k,n); %³õʼÌõ¼þ u(1,:) = exp(x);
%ÇóA¾ØÕó¤Î¶Ô½ÇÔªËØd d = zeros(1,n);
d(1,1) = b1*P1+h*a1; d(2:(n-1),1) = 2*P1; d(n,1) = b2*P1+h*a2;
%ÇóA¾ØÕó¤Î¶Ô½ÇÔªËØÏÂÃæÒ»ÐÐÔªËØe e = -ones(1,n-1); e(1,n-1) = -b2;
%ÇóA¾ØÕó¤Î¶Ô½ÇÔªËØÉÏÃæÒ»ÐÐÔªËØf f = -ones(1,n-1); f(1,1) = -b1;
R = zeros(k,n);%ÇóR %×·¸Ï·¨Çó½â for i = 2:k
R(i,1) = (b1*P2-h*a1)*u(i-1,1)+b1*u(i-1,2)+2*h*c1; for j = 2:n-1
R(i,j) = u(i-1,j-1)+2*P2*u(i-1,j)+u(i-1,j+1); end
R(i,n) = b2*u(i-1,n-1)+( b2*P2-h*a2)*u(i-1,n)+2*h*c2; M = chase(e,d,f,R(i,:)); u(i,:) = M';
plot(x,u(i,:)); axis([0 a0 0 t_max]); pause(0.1) end
output = u;
% »æͼ±È½Ï½âÎö½âºÍÓÐÏÞ²î·Ö½â [X,T] = meshgrid(x,t); Z = exp(X+0.1*T);
surf(X,T,Z),xlabel('x'),ylabel('t'),zlabel('u'),title('½âÎö½â'); figure
Fpg
Fpg
surf(X,T,u),xlabel('x'),ylabel('t'),zlabel('u'),title('ÓÐÏÞ²î·Ö½â');
? ×Ó³ÌÐò2£º
function M = chase(a,b,c,f)
% ×·¸Ï·¨Çó½âÈý¶Ô½Ç¾ØÕ󷽳̣¬Ax=f % aÊǶԽÇÏßϱßÒ»ÐФÎÔªËØ % bÊǶԽÇÏßÔªËØ
% cÊǶԽÇÏßÉϱßÒ»ÐФÎÔªËØ
% MÊÇÇóµÃ¤Î½á¹û£¬ÒÔÁÐÏòÁ¿ÐÎʽ±£´æ n = length(b);
beta = ones(1,n-1); y = ones(1,n); M = ones(n,1);
for i = (n-1):(-1):1 a(i+1) = a(i); end
% ½«a¾ØÕóºÍn¶ÔÓ¦
beta(1) = c(1)/b(1); for i = 2:(n-1)
beta(i) = c(i)/( b(i)-a(i)*beta(i-1) ); end
y(1) = f(1)/b(1); for i = 2:n
y(i) = (f(i)-a(i)*y(i-1))/(b(i)-a(i)*beta(i-1)); end
M(n) = y(n);
for i = (n-1):(-1):1
M(i) = y(i)-beta(i)*M(i+1); end end
? ½á¹û£º
Fpg