|
 
- UID
- 64
- 帖子
- 2885
- 精华
- 13
- 积分
- 6081
- 来自
- 重庆万州
|
洛伦兹方程：
dx/dt=-σ(x-y)，
dy/dt=rx-y-xz，
dz/dt=xy-bz，
的数值求解如何用画板实现，这是画板洛伦兹吸引子可视化的关键。下面是洛仑兹方程数值求解的PASCAL程序，像天书一样，不知谁能解读：
｛Lorenz.PAS｝
Program LorenzAttractorHuajie1993;
uses Graph,Crt,VGAFONT;
var
Gd,Gm，lx,ly,lz,n:integer;
lins1,lins2,sigma,b,r,delta:real;
kx1,kx2,kx3,kx4,ky1,ky2,ky3,ky4,kz1,kz2,kz3,kz4:real;
x,y,z,coef:real;
s:string;
begin
sound(500);delay(200);nosound;TextColor(RED);
x:=2; y:=12; z:=7.03; {initial conditions}
writeln('E.Lorenz Attractor, by Liu Huajie, 1993');
writeln('Input parameter value of "r" ');
writeln('"r" in the scope of ［15,40］, usually it equals 28.');
write(' r=');
readln(r) ;
sigma:=10; b:=8/3; {parameter values}
if registerBGIDriver(@VGADriver) <0 then EXIT;
if RegisterBGIfont(@SansFont) < 0 then EXIT;
if RegisterBGIfont(@TripFont) < 0 then EXIT;
Gd:=Vga; Gm:=Vgahi;Initgraph(Gd,Gm,'');
IF GraphResult <> grOK then HALT(1);
n:=0;delta:=0.0005; {integration step}
settextstyle(1,0,2);SetColor(6);
outtextxy(60,60,'1. Projection of Attractor on Plane YOZ');
outtextxy(60,120,'2. Projection of Attractor on Plane XOZ');
outtextxy(60,180,'3. Projection of Attractor on Plane XOY');
outtextxy(60,260,'Please input number 1,2 or 3,');
outtextxy(60,290,'Then press ENTER key.' );
readln(s);closegraph;
initgraph(gd,gm,'');settextstyle(1,0,3);
outtextxy(10,10,'Press any key to exit!');
rectangle(0,0,639,479);
repeat
begin
n:=n+1;
lins1:=delta*sigma;
kx1:=lins1*(-x+y);
kx2:=kx1-lins1*kx1/2;
kx3:=kx1-lins1*kx2/2;
kx4:=kx1-lins1*kx3;
ky1:=delta*(-x*z+r*x-y);
ky2:=ky1-delta*ky1/2;
ky3:=ky1-delta*ky2/2;
ky4:=ky1-delta*ky3;
lins2:=delta*b;
kz1:=delta*(x*y-b*z);
kz2:=kz1-lins2*kz1/2;
kz3:=kz1-lins2*kz2/2;
kz4:=kz1-lins2*kz3;
x:=x+(kx1+2*kx2+2*kx3+kx4)/6;
y:=y+(ky1+2*ky2+2*ky3+ky4)/6;
z:=z+(kz1+2*kz2+2*kz3+kz4)/6;
lx:=round(x*10);ly:=round(y*8);lz:=round(z*8);
if s='1' then begin
PutPixel(ly+290,460-lz,15); end;
if s='2' then begin
PutPixel(lx+300,460-lz,15); end;
if s='3' then begin
putPixel(lx+290,230-ly,15); end;
end;
until KeyPressed;sound(500);delay(200);nosound;readln;CloseGraph;
end. |
|
