matlab微分方程的求解

function dydt=myode(t,x) m=10; k=100; c=4; dydt=[x(2);-c/m*x(2)-k/m*x(1)]; end [t,y]=ode45(@myode, [0:0.01:20],[1,0]);
1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
3. 利用 Matlab 编程时需要用的主要命令 ode45: 基于显示 Runge-Kutta（4，5）方法求解。对于多数方程来讲，ode45 是最佳的尝试函数命令。 ode23： 基于显示 Runge-Kutta（2，3）方法求解。 ode113： 利用变阶 Adams-Bashforth-Moulton 算法求解。 与 ode45 函数相比， 该方程对于精密度步长及方程难于估计时效更好，但是该方法是多步算法，需要 用到前面几个节点的信息来求解当前节点，效率较低。 4. ode45 调用格式 [t, x]=ode45(@myode, [t0:dt:t1], initial_condition) 从调用格式来看，首先必须生成函数文件“myode” 。 function dydt=myode(t,y) dydt=[fun1; fun2;…;funn]; end 5. 应用实例一：求解微分方程：（1 + y 2 ）y ′ − 2y = 0，初值 y(0)=1 时, t 从 0 至 10 时刻的解。 function dydt=myode(t,y) dydt=[y(1)*2./(1+y(1).^2)]; end [t, y]=ode45(@myode, [0:0.01:10], [1])
0
2
4
6
8
10
12
14
16
18
20
9. 作业：任意找一具有物理背景的二阶或者高阶微分方程或方程组，在给定初 始条件下，编写 matlab 程序求解该微分方程或者方程组。
3
7. 应用实例三： 一质点在空中飞行， 所受的空气阻力方向始终与速度方向相反，大小与速度平方 成正比。求质点的飞行轨迹和飞行路程。 根据题意：我们可以建立如下方程（vx， vy 分别为质点沿 x、y 方向的速度，c 为空气阻力系数） ： dx/dt=vx dy/dt=vy mdvx/dt=-cv^2cosa=-cv*vx mdvy/dt=-cv*vy-mg function dydt=myode(t,x) c=0.01; m=1; g=9.8; dydt=[x(3);x(4);-1/m.*(c.*sqrt((x(3).^2+x(4).^2))).*x(3);-1/m.*(c.*sqrt((x(3). ^2+x(4).^2))).*x(4)-g]; end [t,y]=ode45(@myode,[0:0.01:6],[0;0;50*cos(pi/6);50*sin(pi/6)];
20
10
0
-10
-20
-30
-40
-50
0
20பைடு நூலகம்
40
60
80
100
120
140
2
8. 应用实例四： 如图所示的弹簧－阻尼－质量系统的质量为 m，刚度为 k，阻尼为 c。其自由振 动的微分方程可以表示为：mx + cx + kx = 0，其中 x 为任一时刻质量块离平衡 位置的位移，x 与x分别为位移对时间的一阶、二阶导数，分别称为速度与加速 度。若 m=10, k=100, c=4。编写 Matlab 程序求该质量块在初始位移为 1，初始 速度为 0 时，前 20 秒的位移解。若 c=8,50 时，从新求解上述方程，并画出 c=4,8,50 时的位移曲线图 （提示， 先编写函数文件， 再调用 matlab 命令求解） 。
1
6. 应用实例二： 求解范德蒙方程 y ′′ − 1 − y 2 y ′ + y = 0 的解。时间 0 至 20，初值 y(0)=2, y’(0)=0。 function dydt=myode(t,y) dydt=[y(2)； （1-y(1)^2）*y(2)-y(1)]; end [t, y]=ode45(@myode, [0:0.01:20], [2; 0])
Matlab 程序设计： 微分方程求解 主讲人： 王佐才 1. 引言 Matlab 能够求解的微分方程包括：常微分方程的初值问题，常微分方程的边值 问题，时变常微分方程的初值问题，以及偏微分方程。 2. 常微分方程的初值问题 Matlab 可以求解的常微分方程包括： 显示常微分方程： y ′ = f(t, y) 线性隐式常微分方程：M(t, y)y ′ = f(t, y), 其中，M(t, y) 为矩阵。 全隐式常微分方程： f(t, y, y ′ )=0