|
以下是一段解偏微分方程组的程序,谁能告诉我用的是什么方法解的,我怎么查了好多资料都看不出来啊,谢谢大家了!
d B! {+ m: O' E+ ~: U; C2 h
6 U. X! x o4 q0 \1 |! f! a% % Programme pour leopard % 3 P$ a) m* @: h: I% N/ k+ t
global d K F N
2 o- ~: ^/ a. _N = 64 & }$ g3 p5 _/ Q5 X% y
dx = 1/N; , k' F5 ~ \/ _( [& |0 H+ S
d=10 K=0.01 F=0.1 : W" @3 V" s1 @) i7 a/ u1 |1 N3 ]. |
dt=0.001;
; t/ N5 b% H5 u2 XU=rand(N); V=rand(N);
) @- X6 V9 g" S& d$ sFourU=fft2(U); FourV=fft2(V);
- U) W1 j/ A R% ?3 S4 J( Bfor i=1:N, for j=1:N, i0=i-1; j0=j-1; if i0 > N/2 i0=1+N-i; end if j0 > N/2 j0=1+N-j; end FourU(i,j)=FourU(i,j)/(1+(i0*i0+j0*j0)/8); FourV(i,j)=FourV(i,j)/(1+(i0*i0+j0*j0)/8); end end
+ i9 O# v& n+ X; N7 D, rU=real(ifft2(FourU)); V=real(ifft2(FourV)); 1 f/ `) Y) t3 F8 v: i7 l/ [) a3 a
compteur=0; ' X% f- T9 c3 L" e
for t=0:dt:1, # E- N7 D' a; s' z" a) r
format long max(max(U)) max(max(V)) min(min(V)) 6 t3 o8 o X2 D7 B G' `
FourU=fft2(U); FourV=fft2(V); $ b/ T/ M8 D5 n+ N
for i=1:N, for j=1:N, i0=i-1; j0=j-1; if i0 > N/2 i0=1+N-i; end if j0 > N/2 j0=1+N-j; end FourU(i,j)=FourU(i,j)*exp(-dt*(i0*i0+j0*j0)/10); FourV(i,j)=FourV(i,j)*exp(-d*dt*(i0*i0+j0*j0)/10); end end
/ _9 c. x* P# p8 w0 O1 z. u( fU=real(ifft2(FourU)); V=real(ifft2(FourV)); q f' f5 g5 M
for i=1:N, for j=1:N, P(i,j)=-U(i,j)*V(i,j)^2+F*(1-U(i,j)); G(i,j)=U(i,j)*V(i,j)^2-(F+K)*V(i,j); end end
1 I3 K" N3 h8 L4 ?' b1 L6 uU = U + dt*P; V = V + dt*G; 4 f! u' k* d& E1 I0 G b
for i=1:N, for j=1:N, W(i,j)=V(i,j); M(i,j)=U(i,j); end end ( `- N0 R/ b+ V; I8 M# ?/ J
if rem(compteur,10) == 0 clf; if W>0.6 image(100*W-60); else image(100*W) end if M>0.6 image(100*M-60); else image(100*M); end drawnow; end
+ c% T, W" l' D6 ccompteur = compteur +1; - V7 @7 ^" S0 Y L2 p) F
end %------------------------------------------
|