龙格—库塔法分析Lorenz方程课程设计.doc

数值分析课程设计实验报告龙格库塔法分析Lorenz方程 一、 问题叙述考虑著名的Lorenz方程其中s，r，b为变化区域内有一定限制的实参数，该方程形式简单，表面上看并无惊人之处，但由该方程揭示出的许多现象，促使“混沌”成为数学研究的崭新领域，在实际应用中也产生了巨大的影响。二、 问题分析Lorenz方程实际上是一个四元一阶常微分方程，用解析法精确求解是不可能的，只能用数值计算，最主要的有欧拉法、亚当法和龙格- 库塔法等。为了得到较高精度的，我们采用经典四阶龙格库塔方法求解该问题。三、 实验程序及注释（1）算法程序function T=Runge_Kutta(f,x0,y0,h,n) 定义算法，其中f为待解方程组，x0是初始自变量，y0是初始函数值，h是步长，n为步数if nargin5 n=100; 如果输入参数个数小于5，则步数n=100endr=size(y0);r=r(1); 返回初始输出矩阵的行列数，并将值赋给r(1)s=size(x0);s=s(1); 返回初始输入矩阵的行列数，并将值赋给s(1)r=r+s;T=zeros(r,n+1); T(:,1)=y0;x0; for t=2:n+1 以下是具体的求解过程 k1=feval(f,T(1:r-1,t-1); k2=feval(f,k1*(h/2)+T(1:r-1,t-1);x0+h/2); k3=feval(f,k2*(h/2)+T(1:r-1,t-1);x0+h/2); k4=feval(f,k3*h+T(1:r-1,t-1);x0+h); x0=x0+h; T(:,t)=T(1:r-1,t-1)+(k1+k2*2+k3*2+k4)*(h/6);x0;end（2）主程序function dy=fun(x) 定义函数s=10.0; 给参数s，r，b赋值r=28.0;b=8.0/3;dy(1)=s*(x(2)-x(1); Lorenz方程表达式dy(2)=(r*x(1)-x(3)*x(1)-x(2);dy(3)=x(1)*x(2)-b*x(3);dy=dy;（2）运行程序 T=Runge_Kutta(fun,0,10;10;10,0.01,5000); 调用前面的算法程序plot3(T(1,:),T(2,:),T(3,:); 显示三分量的关系图axis(-20 20 -50 50 0 50) 定义坐标轴长度view(3) 设定观察角度四、 实验数据结果及分析（1）各初始变量相同时的图像分析各初始变量取相同的值10，10，10，运行上述程序后，得到如下图像：从图中可以看出，各初始变量相同时，曲线总是被吸引回奇怪吸引子附近作来回跳跃。初始变量值取为-10，-10，-10 ，20，20，20时，依然如此。图像如下：-10，-10，-10 20，20，20（2）初始值的每个分量变化对图像的影响y分量：0，2，0 0，5，00，15，0 0，20，0从上面可以看出，随着初始y值的增大，奇怪吸引子中曲线在其附近来回跳跃的两个位置中的一个吸引力变弱，另一个吸引力变强。初始y继续增大到某一特定值，情况又会变回来。这说明在空间存在一些区域，当初始位置位于这些区域外时解将出现奇怪吸引子的性质，而在这些区域以内解将呈现普通吸引子的性质。z分量：0，0，20从上图可以看出解的曲线为一直线，这可以从方程的角度来解释。当x=0,y=0时在方程中dx/dt=0,dy/dt=0，x，y 方向的值不发生变化，仅z方向的值变化，因此解为一直线。（3）调整参数r、s、b对图像的影响为便于分析，我们只调整r、s、b三个参数中的任意一个。当只调整b且将初始变量取为0，eps，0。具体情况如下：s=10.0，r=28.0，b=8.0/3 s=10.0，r=28.0，b=9.6/3s=10.0，r=28.0，b=11.0/3 s=10.0，r=28.0，b=15.0/3增大b值时，Lorenz曲线在其附近来回跳跃的两个位置会一个加强，一个减弱。当b达到某一值时，个位置丧失吸引力，另一位置则将曲线完全吸引过来变成普通吸引子。改变s和r的值也有类似的现象。五、 实验结论本实验利用龙格库塔法对Lorenz方程进行了分析，从实验中我们得出，Lorenz方程的解对初始变量和参数r、s、b具有很强的敏感性。数值分析课程设计实验报告数值分析试验报告矩阵的LU分解1. 题目：求4阶矩阵的LU分解2. 方法：杜里特尔分解法3. 程序：function f=LU_decom(A)m,n=size(A)L=eye(n);U=zeros(n);flag=ok;for i=1:n U(1,i)=A(1,i); endfor r=2:n L(r,1)=A(r,1)/U(1,1); end for i=2:n for j=i:n z=0; for r=1:i-1 z=z+L(i,r)*U(r,j); end U(i,j)=A(i,j)-z; end if abs(U(i,i) LU_decom(A)m = 4n = 4L = 1 0 0 0 2 1 0 0 1 2 1 0 3 3 2 1U = 2 4 2 6 0 1 2 3 0 0 3 6 0 0 0 15.拓展在编写程序过程中由于角标较多因此在运行过程中出现了不少角标不对的错误题目；给出函数f(x)=1/(1+25x2)，求f（x在-1,1上取5个和9个等距节点，做最小二乘拟合，得出均方误差方误差。五个节点时，matlab编码为：首先建立M文件，并保存function y=f(x)y=1/(1+25*x2);endx=-1 -0.5 0 0.5 1;for i=1:5y(i)=f(x(i);enda=polyfit(x,y,3)syms xf1=a(1)*x3+a(2)*x2+a(3)*x+a(4) x=-1 -0.5 0 0.5 1;for i=1:5E(i)=(f(x(i)-(a(1)*x(i)3+a(2)*x(i)2+a(3)*x(i)+a(4)2;endsum(E)输出结果为a = -0.0000 -0.6063 -0.0000 0.5737 f1 = -4878849915647781/1298074214633706907132624082305024*x3-1600/2639*x2-5348847520430703/649037107316853453566312041152512*x+1514/2639 （拟合的多项式） ans = 0.3534（均方误差）九个点的时候，matlab编码为：x=-1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1;for i=1:9y(i)=f(x(i);enda=polyfit(x,y,3)syms xf2=a(1)*x3+a(2)*x2+a(3)*x+a(4) x=-1 -0.75 -0.5 -0.25 0 0.25 0.5 0.75 1;for i=1:5E1(i)=(f(x(i)-(a(1)*x(i)3+a(2)*x(i)2+a(3)*x(i)+a(4)2;endsum(E1)输出结果为：a = -0.0000 -0.5609 0.0000 0.4855 f2 = -728732707776039/2535301200456458802993406410752*x3-1263030908712921/2251799813685248*x2+4915246442354361/20282409603651670423947251286016*x+1093229300764671/2251799813685248 （最小二乘拟合多项式）ans =0.3350（均方误差） 用复合梯形公式求积分的值。function i=combinetraprl(f,a,b,eps)%复化梯形公式求函数f在区间a，b上的定积分%函数名：f%积分下限：a%积分上限：b%积分精度：eps%积分值：i%积分划分的子区间个数：stepn=1;h=(b-a)/2;i1=0;i2=(subs(sym(f),findsym(sym(f),a)+subs(sym(f),findsym(sym(f),b)/h;while abs(i2-i1)eps n=n+1; h=(b-a)/n; i1=i2; i2=0; for i=0:n-1 x=a+h*i; x1=x+h; i2=i2+(h/2)*(subs(sym(f),findsym(sym(f),x)+. subs(sym(f),findsym(sym(f),x1); endend四阶龙格-库塔法分别求解下列初值问题；functionx,y=runge_kutta(fun,x0,xt,y0,pointnum,varargin)if nargin3 y0=0;endy(1,:)=y0(:);h=(xt-x0)/(pointnum-1);x=x0+0:pointnum*h;for k=1:pointnum f1=h*feval(fun,x(k),y(k,:),varargin:); f1=f1(:); f2=h*feval(fun,x(k)+h/2,y(k,:)+f1/2,varargin:); f2=f2(:); f3=h*feval(fun,x(k)+h/2,y(k,:)+f2/2,varargin:); f3=f3(:); f4=h*feval(fun,x(k)+h,y(k,:)+f3,varargin:); f4=f4(:); y(k+1,:)=y(k,:)+(f1+2*(f2+f3)+f4)/6;endfunctionx,y=euler(fun,x0,xt,y0,pointnum)if narginx0=0;xt=1;fun=inline(9*y/(1+2*x),x,y);y0=1;pointnum=10;x1,yr=runge_kutta(fun,x0,xt,y0,pointnum)一 运行结果：1） Runge_kutta方法：x1 =0 0.1111 0.2222 0.3333 0.4444 0.5556 0.6667 0.7778 0.8889 1.0000 1.1111yr =1.0000 2.4611 5.2136 9.9212 17.4198 28.7302 45.0700 67.8650 98.7600 139.6286 192.5834小行星轨道问题：一天文学家要确定一颗小行星绕太阳运行的轨道，他在轨道平面内建立以太阳为原点的直角坐标系，在五个不同的对小行星作了五次观察，测得轨道上五个点的坐标数据（单位：万公里）如下表所示：P1P2P3P4P5X坐标5360558460628596666268894Y坐标602611179169542349268894由开普勒第一定律知，小行星轨道为一椭圆，椭圆的一般方程可表示为：a1x2+2a2xy+a3y2+2a4x+2a5y+1=0现需要建立椭圆的方程以供研究。（1）分别将五个点的数据代入椭圆一般方程中，写出五个待定系数满足的等式，整理后写出线性方程组AX=R以及方程组的系数矩阵和右端项b；（2）用MARLAB求低阶方程的指令Ab求出待定系数a1,a2,a3,a4,a5.答：1)2)一源程序：A=53