matlab 洛伦兹方程,求解洛伦兹方程的MATLAB程序

求解洛伦兹方程的MATLAB程序

clear;a=0;b=100;h1=0.0001;h2=0.1;ya=5;2;10;sigma=10;gamma=28;rho=8/3;z=f2(sigma,gamma,rho);y1,x1=Euler(z,a,b,h1,ya);y2,x2=Rungkuta4(z,a,b,h2,ya);y3,x3=Rungkuta4(z,a,b,h2/2,ya);x1=x1x2=x2y1=y1;y2=y2y3=y3z1=sqrt(sum(y2(end,:)-y1(end,:).2)z2=sqrt(sum(y3(end,:)-y1(end,:).2)figureplot(x,y(:,1)title(y1随时间t的变化情况);xlabel(t);ylabel(数值解y);figureplot(x,y(:,2)title(y2随时间t的变化情况);xlabel(t);ylabel(数值解y);figureplot(x,y(:,3)title(y3随时间t的变化情况);xlabel(t);ylabel(数值解y);figureplot3(y(:,1),y(:,2),y(:,3);title(y1,y2,y3的三维运动轨道);xlabel(y1);ylabel(y2);zlabel(y3)function y x=Euler(f,a,b,h1,ya)x=a:h1:b;N=length(x);y=zeros(3,N);y(:,1)=ya;for i=1:N-1y(:,i+1)=y(:,i)+h1*feval(f,x(i),y(:,i);endfunction y,x=Rungkuta4(f,a,b,h2,ya)x=a:h2:b;N=length(x);y=zeros(3,N);y(:,1)=ya;for i=1:N-1k1=feval(f,y(:,i);k2=feval(f,y(:,i)+(h2/2)*k1);k3=feval(f,y(:,i)+(h2/2)*k2);k4=feval(f,y(:,i)+h2*k3);y(:,i+1)=y(:,i)+(h2/6)*(k1+2*k2+2*k3+k4);end