大学数学实验四_常微分方程数值解

输出结果为：
MATLAB 自带龙格-库塔公式： n=1000; h=60/n; t=0:h:60; v(1)=0; [t,v]=ode45(@rocket,t,v(1)); z(1)=0; for i = 1:n
a = z(i); b = v(i); c = v(i+1); a = a + (b+c)/2*h; z(i+1) = a; end; [t(n+1),z(n+1)] hold on; grid on; plot (t,z,'r') %输出时间、高度
令
则
故有
10
故有
两边积分得
故
取倒数得
故有
即
将 代入得
至此，小船渡河路径的解析式已求得。 用 MATLAB 求解小船渡河的路径： function dy=duhe(x,y) v1=1; v2=2; k=v1/v2; dy=k*((y/(100-x)).^2+1).^0.5-(y/(100-x)); x=0:1:100; y(1)=0; [x,y]=ode23(@duhe,x,y(1)); hold on; grid on;
前 60s 火箭的 v – t 图如下：
t = 60s 时的火箭的速度为（经典龙格-库塔公式）：
n = 10
n = 50
n = 100
265.0392
268.0522
267.4982
n = 5000
n = 10000
n = 50000
266.9538
266.9483
266.9438
使用 MATLAB 自带的龙格-库塔方法计算： function dv=rocket(t,v) dv=(32000-0.4*v.^2)/(1400-18*t)-10;
n=100; h=20/n; t=60:h:80; v(1)=266.93; [t,v]=ode45(@rockett,t,v(1)); hold on; grid on; plot(t,v,'r');
7
得到结果如下：
火箭的速度-时间曲线（单位：米/秒，秒）
下面求高度随时间变化的曲线。 之前已经得到前 60 秒的图线，现在求 60 秒之后的图线。 n=1000; h=20/n; t=60:h:80; v(1)=266.93; [t,v]=ode45(@rockett,t,v(1)); z(1)=12131.48; f(1)=0.934; for i = 1:n
前 60s 火箭的 v – t 图如下：
t = 60s 时的火箭的速度为（向前欧拉公式）：
n = 10
n = 50
n = 100
273.2605
268.0559
267.5001
n = 5000
n = 10000
பைடு நூலகம்
n = 50000
266.9539
266.9483
266.9438
n = 500 267.0543 n = 100000 266.9432
本题建模时有如下几个假设： （1）燃料燃烧率为一恒定值； （2）空气阻力的大小严格正比于速度的平方； （3）火箭能够达到的最高高度不是太高，故重力加速度可视为恒定值 g = 10 kg·m/s2； （4）燃料燃烧产生的推力恒定； （5）引擎关闭时没有延迟时间。 本章学习了用欧拉方法和龙格-库塔方法求解常微分方程。MATLAB 提供了用龙格-库塔方法 求解的语句 ode45(@f, x, y0)。求解引擎关闭瞬间的速度、高度、加速度时，我分别用了向 前欧拉公式、经典的龙格-库塔公式和 ode45 语句 3 种不同的方法求解，以比较所得结果的 差异。之后物理过程的求解都只使用了 MATLAB 自带的 ode45 语句。 燃料耗尽时，引擎关闭，此时 t = 60s。先讨论 60s 以内的情况。 根据牛顿第二定律，F 推 – f = m(a + g)，故 32000 – 0.4v2 = (1400 – 18t)(a +g)。 根据上述分析，有
1
根据向前欧拉公式编程解第一小问： n=10; h=60/n; %步长 t=0:h:60; v=0:h:60; %设定v(i)的初值,v(1)=0,其余的v(i)值在后面都会另予赋值 for i=1:n
a = v(i); a = a +((32000-0.4*a.^2)/(1400-18*h*i)-10)*h; %向前欧拉公式 v(i+1) = a; end; hold on; grid on; plot(t,v,'r'); %输出v-t图线 [t(n+1)',v(n+1)'] %输出时间、速度
n = 500 267.0539 n = 100000 266.9432
n = 1000 266.9983 n = 1000000 266.9427
3
n=1000; h=60/n; t=0:h:60; v(1)=0; [t,v]=ode45(@rocket,t,v(1)); a(1)=90/7; for i = 1:(n+1)
a = v(i); b = a +((32000-0.4*a.^2)/(1400-18*h*i)-10)*h; v(i+1) = b; c =c + (a+b)*h/2; e = t(i+1); d =(32000-0.4*b.^2)/(1400-18*e)-10; z(i+1) = c; f(i+1) = d; end; hold on; grid on; plot(t,z,'r'); [t(n+1)',z(n+1)', f(n+1)'] %输出时间、高度、加速度
n = 1000 266.9985 n = 1000000 266.9427
编写用经典龙格-库塔公式解第一小问： n=10; h=60/n; %步长 t=0:h:60; v=0:h:60; %设定v(i)的初值,v(1)=0,其余的v(i)值在后面都会另予赋值 for i=1:n
a = v(i);
2
5
d = a + h*k3; k4 = (32000-0.4*d.^2)/(1400-18*x1)-10; e = a +(k1+2*k2+2*k3+k4)*h/6; v(i+1) = e; p = t(i+1); q =(32000-0.4*e.^2)/(1400-18*p)-10; g(i+1) = q; s = s + (a+e)*h/2; z(i+1) = s; end; [t(n+1)',z(n+1)', g(n+1)'] %输出时间、高度、加速度
大学数学实验四 常微分方程数值解 实验报告
【实验目的】 1、掌握用 MATLAB 软件求微分方程初值问题数值解的方法。 2、通过实例学习用微分方程模型解决简化的实际问题。 3、 了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。
【实验内容】
2 小型火箭初始质量为 1400kg，其中包括 1080kg 燃料。火箭竖直向上发射时，燃料燃烧率 为 18kg/s，由此产生 32000N 推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正 比于速度的平方，比例系数为 0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火 箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
6
输出结果为
故，当 t =60s 时，火箭的速度、高度、加速度为
v (m/s)
向前欧拉公式
266.94
经典龙格-库塔公式
266.94
ode45(@rocket,t,v)
266.93
h (m) 12131 12131 12131
从上表可以看出，用三种不同方法求得的结果相差不大。
60s 后，火箭作竖直上抛运动，加速度为
a (m/s2) 0.9266 0.9266 0.9343
下面用 MATLAB 求火箭速度、高度、加速度随时间变化曲线，并求火箭到达的最高高度。 先求速度随时间变化的曲线。前 60s 的曲线已经得到，下面做 60 秒之后的曲线。 function dv=rockett(t,v) dv=-0.4*v.^2/320-10;
4
输出的结果为
输出的火箭的 h – t 图线如下：
经典龙格-库塔公式： n=100000; h=60/n; t=0:h:60; v(1)=0; z(1)=0; g(1)=90/7; for i=1:n
a = v(i); s = z(i); x1 = h*i; k1 = (32000-0.4*a.^2)/(1400-18*x1)-10; x1 = x1 + h/2; b = a + h/2*k1; k2 = (32000-0.4*b.^2)/(1400-18*x1)-10; c = a + h/2*k2; k3 = (32000-0.4*c.^2)/(1400-18*x1)-10; x1 = x1 + h/2;
c = v(i+1); e = t(i+1); d =(32000-0.4*c.^2)/(1400-18*e)-10; f(i+1) = d; end; hold on; grid on; plot(t,f,'r');
9
再作60秒之后的加速度曲线： 作高度-时间曲线时已有说明（第8页）。
火箭的加速度-时间曲线（单位：米/秒 2，秒） 实验结果： （1）火箭的引擎关闭的瞬间，t =60s，v =266.9m/s，h =12131m，a =0.93m/s2。 （2）火箭在 71 秒时到达最高高度 13050 米，此时加速度为– 10 m/s2（即重力加速度）。 5 一只小船渡过宽为 d 的河流，目标是起点 A 正对着的对岸的 B 店。已知河水流速为 v1， 船在静水中的速度为 v2，且 v1:v2=k。（1）建立描述小船航线的数学模型并求其解析解；（2） 设 d=100m，v1=1m/s，v2=2m/s，用数值解法求渡河所需时间、任意时刻小船的位置和航行 曲线，作图并和解析解作比较；（3）若流速 v1=0, 0.5, 1.5, 2m/s，结果将如何？ 建立数学模型并求解析式： 以 A 为原点，v1 方向为 y 轴正方向，B 点的坐标为(100, 0)。 【由于第三小题中 d 的值并未改变，故在此默认 d=100m。事实上，如果 d 取其他值，只需 将下面求解得到的解析式中的 100 换成 d 即可。】 设小船由 A 到 B 的轨迹为函数 y = f(x)，那么有
x1 = h*i; k1 = (32000-0.4*a.^2)/(1400-18*x1)-10; x1 = x1 + h/2; b = a + h/2*k1; k2 = (32000-0.4*b.^2)/(1400-18*x1)-10; c = a + h/2*k2; k3 = (32000-0.4*c.^2)/(1400-18*x1)-10; x1 = x1 + h/2; d = a + h*k3; k4 = (32000-0.4*d.^2)/(1400-18*x1)-10; a = a +(k1+2*k2+2*k3+k4)*h/6; %经典龙格-库塔公式 v(i+1) = a; end; hold on; grid on; plot(t,v,'k'); %输出v-t图线 [t(n+1)',v(n+1)'] %输出时间、速度
b = v(i); e = t(i); d =(32000-0.4*b.^2)/(1400-18*e)-10; a(i) = d; end; [t(n+1),v(n+1),a(n+1)] %输出时间、速度、加速度
输出结果为
下面求引擎熄灭时火箭的高度和加速度。 向前欧拉公式： n=100000; h=60/n; t=0:h:60; v=0:h:60; c = 0; z(1)=0; f(1)=90/7; for i=1:n
8
输出的曲线如下。
火箭的高度-时间曲线（单位：米，秒） 输出的数据如下（部分）：
故火箭在 71 秒时到达最高高度 13050 米，此时加速度为– 10 m/s2。 下面求加速度随时间变化的曲线。先作前 60 秒的加速度曲线： n=1000; h=60/n; t=0:h:60; v(1)=0; [t,v]=ode45(@rocket,t,v(1)); f(1)=90/7; for i = 1:n
a = z(i); b = v(i); c = v(i+1); a = a + (b+c)/2*h; z(i+1) = a; d =-0.4*b.^2/320-10; f(i+1) = d; end; hold on; grid on; plot(t,z,'r'); %后面画加速度曲线时，将本句中的z改为f即可 [t,z',v,f'] %输出时间、高度、速度、加速度