洛伦兹方程是一组典型的非线性方程，来源于大气物理研究。在一定的参数下，方程的解表现出了非常复杂而又有趣的现象。本专题要求采用数值积分的方法，求解如下洛伦兹方程： 其中，。 为从0 开始的时间变量，积分初值为任意的不全为零的数值即可。 要求：（1）编写程序求解上述方程（采用C/C++、Matlab/SImulink 均可）；（2）绘制三维运动轨道；（3）观察不同初值的影响，观察最终运动状态。(分组进行，同以前的专题研究要求一样) 提交：（1）程序源码；（2）1~2 幅三维截图；（3）描述你的心得和体会。 提示：可以采用欧拉法或 4 阶龙格-库塔法等数值积分的方法。 我们用MATLAB 编程，基于4 阶 Runge-Kutta 方法的迭代作出了洛伦玆吸引子的求解 程序主体： a=1:1:3000;b=1:1:3000;c=1:1:3000; x1=20;y1=30;z1=40; x2=0;y2=0;z2=0; h=0.1; kfx=0;ksx=0;ktx=0;kfhx=0;kfy=0;ksy=0;kty=0;kfhy=0;kfz=0;ksz=0;ktz=0;kfhz=0; for i=1:1:3000 abx=sym('10*(y0-x0)'); aby=sym('x0*(28-z0)-y0'); abz=sym('x0*y0-8/3*z0'); x0=x1;y0=y1;z0=z1; kfx=eval_r(abx);kfy=eval_r(aby);kfz=eval_r(abz); x0=x1+1/2*h*kfx;y0=y1+1/2*h*kfy;z0=z1+1/2*h*kfz; ksx=eval_r(abx);ksy=eval_r(aby);ksz=eval_r(abz); x0=x1+1/2*h*ksx;y0=y1+1/2*h*ksy;z0=z1+1/2*h*ksz; ktx=eval_r(abx);kty=eval_r(aby);ktz=eval_r(abz); x0=x1+1/2*h*ktx;y0=y1+1/2*h*kty;z0=z1+1/2*h*ktz; kfhx=eval_r(abx);kfhy=eval_r(aby);kfhz=eval_r(abz); x2=x1+h/6*(kfx+2*ksx+2*ktx+kfhx); y2=y1+h/6*(kfy+2*ksy+2*kty+kfhy); z2=z1+h/6*(kfz+2*ksz+2*ktz+kfhz); a(i)=x2;b(i)=y2;c(i)=z2; x1=x2;y1=y2;z1=z2; end plot3(a,b,c) 截图为 这个是我们用MatLab 的数值解法解微分方程组的方程得出来的结果。这个源程序是 fu nction y =lu c(t,x ) y =[10*(x (2)-x (1));x (1)*(28-x (3))-x (2);x (1)*x (2)-8/3*x (3)]; end %这个为所写的函数 [t,x ]=ode45('lu c',[0,300],[1,2,3]); % 调用函数 plot3(x (:,1),x (:,2),x (:,3)) 心得体会：通过这个研究我们学会了如何查找资料，如何去选择资料，以及如何去利用资料来描述我们的问题，并试图解决它。我们在这里遇到了一些难题，比如如何将数学理论转化为程序，我们是先把数学思想吃透，然后就我们将其数学思想转化为程序的各个流程，然后在电脑上调试，由于这个程序的ho ld o n 比较多，所以运行起来速度比较慢，我们只好耐心的等待。通过这次专题研究我们几个也对团体的合作分工有了更进一步的了解，就是查资料的时候分工查，你查网络资料，我查图书资料等。然后再一起来汇总来分析哪些资料有用，哪些资料是关键资料。这样我们配合的很好。
