2020年计算物理解常微分方程编程 .pdfVIP

例 1： function jyjs [T,Y]=ode45(yjs,[0:10],[2,1]); plot(T,Y(:,1),T,Y(:,2)) function ydot=yjs(t,y) ydot=[y(2); 4 ]; 例 2：研究受空气阻尼的抛体运动 function znxp global m d e % 用全局变量传递参数 m=1;b=[0,0.2,0.2];p=[0,0,1]; % 传递给方程的参数 px=[4.6;4.5;4.5]; % 说明文字的 x 坐标 py=[3.5;1.8;0.4]; % 说明文字的 y 坐标 strdd{1}= 无阻尼 ; % 图形说明文字 strdd{2}= 阻力正比于 v; strdd{3}= 阻力正比于 v^2; figure for i=1:3 % 三次循环计算三种条件下的运动 d=b(i); e=p(i); % 选定传递参数的值 [t,y]=ode45(znxpfun,[0:0.01:10],[0,3,0,5]); H{i}=max(y(:,3)) % 轨迹最高点 T{i}=t(find(y(:,3)==H{i})) % 到达最高点时间 vx0{i}=y(find(y(:,3)==H{i}),2) % 到最高点速度 subplot(2,1,1) % 分区画图 axis([0 6 -70 2]); % 确定坐标围 hold on % 保留图形 xlabel(x); ylabel(y); %x,y 轴的标注 plot(y(:,1),y(:,3)); % 粒子的空间轨迹 subplot(2,1,2) axis([0 10 0 4]) hold on xlabel(t); ylabel(dx/dt) text(px(i),py(i),strdd{i}); % 加注解文字 plot(t,y(:,2)) % 速度随时间变化图 end % function ydot=znxpfun(t,y) % 要解的常微分方程 global m d e ydot=[y(2); -d/m*y(2)*(y(2).^2+y(4).^2)^(e/2); y(4); -9.8-d/m*y(4)*(y(2).^2+y(4).^2)^(e/2)]; 例 3: 混沌理论的洛仑兹方程 function lor% 建立函数文件 axis([10 50 -50 50 -50 50]) % 设定坐标轴围 view(3) % 设定视角 hold on % 保持设置 title(Lorenz Attractor) % 设定标题 options =odeset(OutputFcn,odephas3) ; % 设定选项 ,绘出解向量中前三个元素序列构成的三 维相空间图 [t,y]=ode23(lorfun,[0,20],[0,0,eps],options); % 解方程 sol=ode23(lorfun,[0,20],[0,0,eps]) % 以结构数组输出解 ss=deval(sol,10) % 查看 t=10 时的解 function ydot=lorfun(t,y) ydot=[-8/3*y(1)+y(2)*y(3); -10*y(2)+10*y(3); -y(2)*y(1)+35*y(2)-y(3)]; 例 4：小球弹跳 ballode 例 5：求解常微分方程 F=(t,y)[y(2);-y(1)]; %F=(t,y)[0,1;-1,0]*y; %使用矩阵形式来编写匿名函数 %ode45(F,[0,10],[2,0]) sol=ode45(F,[0,10],[0,1]) ss=deval(sol,5.15) options=odeset(outputfcn,odephas2);% 改变原来输出的设置 ode45(F,[0,10],[0,1],options) 例 5：火焰问题