大学数学实验四常微分方程数值解
实验报告
【实验目的】
1、掌握用MATLAB软件求微分方程初值问题数值解的方法。
2、通过实例学习用微分方程模型解决简化的实际问题。
3、了解欧拉方法和龙格-库塔方法的基本思想和计算公式，及稳定性等概念。【实验内容】
2 小型火箭初始质量为1400kg，其中包括1080kg燃料。火箭竖直向上发射时，燃料燃烧率为18kg/s，由此产生32000N推力，火箭引擎在燃料用尽时关闭。设火箭上升时空气阻力正比于速度的平方，比例系数为0.4kg/m，求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点时的高度和加速度，并画出高度、速度、加速度随时间变化的图形。
本题建模时有如下几个假设：
（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)。
根据上述分析，有
根据向前欧拉公式编程解第一小问：
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)'] %输出时间、速度
前60s火箭的v–t图如下：
编写用经典龙格-库塔公式解第一小问：
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);
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)'] %输出时间、速度
前60s火箭的v–t图如下：
使用MATLAB自带的龙格-库塔方法计算：
function dv=rocket(t,v)
dv=(32000-0.4*v.^2)/(1400-18*t)-10;
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)
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
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)'] %输出时间、高度、加速度
输出的结果为
输出的火箭的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;
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)'] %输出时间、高度、加速度输出结果为：
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') %输出时间、高度
输出结果为
故，当t =60s时，火箭的速度、高度、加速度为
从上表可以看出，用三种不同方法求得的结果相差不大。
60s 后，火箭作竖直上抛运动，加速度为
下面用MATLAB求火箭速度、高度、加速度随时间变化曲线，并求火箭到达的最高高度。先求速度随时间变化的曲线。前60s的曲线已经得到，下面做60秒之后的曲线。
function dv=rockett(t,v)
dv=-0.4*v.^2/320-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');
得到结果如下：
火箭的速度-时间曲线（单位：米/秒，秒）
下面求高度随时间变化的曲线。
之前已经得到前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
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'] %输出时间、高度、速度、加速度
输出的曲线如下。
火箭的高度-时间曲线（单位：米，秒）
输出的数据如下（部分）：
故火箭在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
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');
再作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)，那么有
令
则
故有
故有
两边积分得
故
取倒数得
故有
即
将
代入得
至此，小船渡河路径的解析式已求得。
用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;
plot(x,y,'b');
[x,y]
输出的图线如下。
小船航行的路径输出的数据如下。
x y
0 0
1.0000 0.4975
2.0000 0.9900
3.0000 1.4773
4.0000 1.9596
5.0000 2.4367
6.0000 2.9087
7.0000 3.3754
8.0000 3.8369
9.0000 4.2930
10.0000 4.7437
11.0000 5.1890
12.0000 5.6288
13.0000 6.0631
14.0000 6.4919
15.0000 6.9150
16.0000 7.3326
17.0000 7.7444
18.0000 8.1504
19.0000 8.5506
20.0000 8.9450 21.0000 9.3333
22.0000 9.7157
23.0000 10.0920
24.0000 10.4621
25.0000 10.8261
26.0000 11.1839
27.0000 11.5354
28.0000 11.8805
29.0000 12.2191
30.0000 12.5511
31.0000 12.8765
32.0000 13.1952
33.0000 13.5071
34.0000 13.8122
35.0000 14.1103
36.0000 14.4015
37.0000 14.6855
38.0000 14.9624
39.0000 15.2319
40.0000 15.4939
41.0000 15.7484
42.0000
15.9951
43.0000 16.2342
44.0000 16.4653
45.0000 16.6886
46.0000 16.9038
47.0000 17.1108
48.0000 17.3094
49.0000 17.4994
50.0000 17.6807
51.0000 17.8531
52.0000 18.0165
53.0000 18.1706
54.0000 18.3155
55.0000 18.4510
56.0000 18.5768
57.0000 18.6928
58.0000 18.7986
59.0000 18.8939
60.0000 18.9786
61.0000 19.0523
62.0000 19.1147
63.0000 19.1657
64.0000 19.2051
65.0000 19.2329
66.0000 19.2486
67.0000 19.2515
68.0000 19.2413
69.0000 19.2173
70.0000 19.1792
71.0000 19.1265
72.0000 19.0585
73.0000 18.9749
74.0000 18.8754
75.0000 18.7595
76.0000 18.6263
77.0000 18.4746
78.0000 18.3035
79.0000 18.1120
80.0000 17.8991
81.0000 17.6641
82.0000 17.4060
83.0000 17.1227
84.0000 16.8121
85.0000 16.4721
86.0000 16.1006
87.0000 15.6960
88.0000 15.2544
89.0000 14.7713
90.0000 14.2420
91.0000 13.6621
92.0000 13.0231
93.0000 12.3143
94.0000 11.5246
95.0000 10.6325
96.0000 9.6113
97.0000 8.4111
98.0000 6.9394
99.0000 4.9582 100.0000 0.0000
下面比较解析解和数值解。
x=0:1:100;
y(1)=0;
[x,y]=ode23(@duhe,x,y(1)); z(1)=0;
for i=2:101
a = x(i);
b = 100-a;
c = 1/(100.^0.5);
d = b.^1.5;
e = b.^(-0.5);
f = 1/(2*c*e)-c*d/2;
z(i)=f; %解析解
end
hold on;
grid on;
plot(x,y,'b')
plot(x,z,'r--');
输出结果如下。
可见，数值解和解析解之间十分吻合。为了了解它们之间的误差究竟有多大，求它们之间的差值对x轴的曲线和相对误差对x轴的曲线。
clear all;
n=10000;
h=100/n;
x=0:h:100;
y(1)=0;
[x,y]=ode23(@duhe,x,y(1));
z(1)=0;
z2(1)=0;
for i=2:(n+1)
a = x(i);
m = y(i);
b = 100-a;
c = 1/(100.^0.5);
d = b.^1.5;
e = b.^(-0.5);
f = 1/(2*c*e)-c*d/2;%精确值
g = m-f;%绝对误差
s = g/f;%相对误差
z(i) = g;
z2(i) = s;
end
subplot(1,2,1);
hold on;
grid on;
plot(x,z,'r');
subplot(1,2,2);
hold on;
grid on;
plot(x,z2,'b');
输出的曲线如下。
从上图中观察到，数值解与解析解之间的绝对误差很小，而相对误差随着x的增大有所增大（这可能是由于局部截断误差的累积造成的），但数值一直都很小，最大也没有超过0.35%，可认为数值解是可靠的。
用MATLAB求小船渡河的时间：
function dx=duhe2(t,x)
v1=1;
v2=2;
r=sqrt((x(2)).^2+(100-x(1)).^2);
dx=[v2*(100-x(1))/r;v1-v2*x(2)/r];
ts=0:1:70;
x0=[0,0];
[t,x]=ode23(@duhe2,ts,x0);
hold on;
grid on;
plot(t,x); [t,x]
输出的数据如下。
t (s) x (m) y (m)
0 0 0
1.0000
2.0000 0.9899
2.0000
3.9997 1.9595
3.0000 5.9991 2.9082
4.0000 7.9977 3.8357
5.0000 9.9956 4.7415
6.0000 11.9922 5.6254
7.0000 13.9874 6.4867
8.0000 15.9809 7.3249
9.0000 17.9723 8.1395
10.0000 19.9613 8.9299
11.0000 21.9476 9.6960
12.0000 23.9309 10.4373
13.0000 25.9108 11.1531
14.0000 27.8866 11.8428
15.0000 29.8580 12.5058
16.0000 31.8245 13.1414
17.0000 33.7856 13.7490
18.0000 35.7410 14.3284
19.0000 37.6902 14.8793
20.0000 39.6323 15.4007
21.0000 41.5667 15.8921
22.0000 43.4926 16.3526
23.0000 45.4092 16.7814
24.0000 47.3159 17.1778
25.0000 49.2124 17.5419
26.0000 51.0977 17.8734
27.0000 52.9708 18.1713
28.0000 54.8304 18.4347
29.0000 56.6755 18.6628
30.0000 58.5049 18.8546
31.0000 60.3174 19.0092
32.0000 62.1130 19.1275
33.0000 63.8904 19.2093
34.0000 65.6479 19.2535 35.0000 67.3837 19.2592
36.0000 69.0960 19.2253
37.0000 70.7828 19.1510
38.0000 72.4425 19.0351
39.0000 74.0759 18.8799
40.0000 75.6819 18.6859
41.0000 77.2574 18.4520
42.0000 78.7994 18.1771
43.0000 80.3049 17.8601
44.0000 81.7710 17.4999
45.0000 83.1949 17.0955
46.0000 84.5820 16.6526
47.0000 85.9318 16.1732
48.0000 87.2386 15.6557
49.0000 88.4967 15.0988
50.0000 89.7003 14.5009
51.0000 90.8439 13.8606
52.0000 91.9220 13.1764
53.0000 92.9449 12.4568
54.0000 93.9132 11.7055
55.0000 94.8186 10.9216
56.0000 95.6529 10.1041
57.0000 96.4183 9.2554
58.0000 97.1417 8.3870
59.0000 97.8033 7.4968
60.0000 98.3784 6.5814
61.0000 98.8423 5.6372
62.0000 99.2207 4.6729
63.0000 99.5280 3.6968
64.0000 99.7909 2.7113
65.0000 100.0040 1.7217
66.0000 100.0075 0.7234
67.0000 100.0000 0.0000
68.0000 100.0000 0.0000
69.0000 100.0000 0.0000
70.0000 100.0000 0.0000
根据上面这些数据，我们可以知道任意时刻小船在坐标系中的确定位置，并且我们可以知道，当河水流速为1m/s和船速为2m/s时，小船渡河的具体时间为67s。
输出的图形如下。
小船渡河路径上的点与时间的关系曲线
若流速v1=0, 0.5, 1.5, 2m/s，相应的k = 0, 0.25, 0.75, 1。
用MA TLAB输出小船前进的路径如下图所示。
可见，随着k的增大，小船航线的曲率越来越大。
当k =1时，小船不能到达B点。
这一点从解析式中也很好理解。小船如果能到达B点，那么x在100附近时，y对x的导数应该小于0。但是如果k =1或k >1式，根据
可知，y对x的导数恒大于0，故小船不可能到达B点。
v1=0时，小船的位置(x, y)随时间变化的函数x(t)和y(t)的图像如下。此时小船沿着x轴以2m/s 的速度匀速向B点前进。
部分时刻小船的坐标如下。
小船渡河时间为50秒。
v1=0.5时，小船的位置(x, y)随时间变化的函数x(t)和y(t)的图像如下。
任意时刻小船的坐标如下。
t (s) x (m) y (m)
0 0 0
1.0000
2.0000 0.4950
2.0000
3.9999 0.9797
3.0000 5.9998 1.4541
4.0000 7.9994 1.9178
5.0000 9.9989 2.3707
6.0000 11.9981 2.8126
7.0000 13.9969 3.2431
8.0000 15.9952 3.6620
9.0000 17.9930 4.0690
10.0000 19.9903 4.4638
11.0000 21.9868 4.8462
12.0000 23.9826 5.2161
13.0000 25.9775 5.5729
14.0000 27.9714 5.9163
15.0000 29.9642 6.2459
16.0000 31.9556 6.5613
17.0000 33.9457 6.8621
18.0000 35.9342 7.1482
19.0000 37.9212 7.4192
20.0000 39.9062 7.6745
21.0000 41.8891 7.9135
22.0000 43.8697 8.1357
23.0000 45.8478 8.3404
24.0000 47.8231 8.5270
25.0000 49.7956 8.6957
26.0000 51.7649 8.8458
27.0000 53.7306 8.9764
28.0000 55.6923 9.0864 29.0000 57.6495 9.1749
30.0000 59.6019 9.2409
31.0000 61.5488 9.2835
32.0000 63.4904 9.3030
33.0000 65.4260 9.2984
34.0000 67.3546 9.2679
35.0000 69.2752 9.2098
36.0000 71.1870 9.1224
37.0000 73.0889 9.0039
38.0000 74.9804 8.8539
39.0000 76.8604 8.6710
40.0000 78.7268 8.4524
41.0000 80.5779 8.1951
42.0000 82.4129 7.8989
43.0000 84.2295 7.5612
44.0000 86.0231 7.1762
45.0000 87.7894 6.7383
46.0000 89.5254 6.2449
47.0000 91.2251 5.6910
48.0000 92.8807 5.0689
49.0000 94.4822 4.3711
50.0000 96.0136 3.5848
51.0000 97.4502 2.6939
52.0000 98.7475 1.6730
53.0000 99.7914 0.4715
54.0000 100.0000 0.0000
55.0000 100.0000 0.0000
56.0000 100.0000 0.0000
小船渡河时间为54秒。
v1=1.5时，小船的位置(x, y)随时间变化的函数x(t)和y(t)的图像如下。小船的坐标随时间的变化如下。
t (s) x (m) y (m) 0 0 0 2.0000 3.9994 2.9392 4.0000 7.9949 5.7536 6.0000 11.9825 8.4387 8.0000 15.9572 10.9895 10.0000 19.9134 13.4008 12.0000 23.8462 15.6688 14.0000 27.7489 17.7898 16.0000 31.6137 19.7588 18.0000 35.4330 21.5711 20.0000 39.2008 23.2255 22.0000 42.9088 24.7213 24.0000 46.5466 26.0559 26.0000 50.1036 27.2265 28.0000 53.5699 28.2309 30.0000 56.9415 29.0761 32.0000 60.2094 29.7660 34.0000 63.3621 30.3022 36.0000 66.3884 30.6863 38.0000 69.2781 30.9209 40.0000 72.0367 31.0223 42.0000 74.6598 30.9994 44.0000 77.1381 30.8573 46.0000 79.4625 30.6008 48.0000 81.6237 30.2351 50.0000 83.6366 29.7775 52.0000 85.5143 29.2435 54.0000 87.2503 28.6379
56.0000 88.8384 27.9654 58.0000 90.2722 27.2310 60.0000 91.5613 26.4441 62.0000 92.7449 25.6218 64.0000 93.8193 24.7673 66.0000 94.7780 23.8831 68.0000 95.6149 22.9715 70.0000 96.3298 22.0358 72.0000 96.9908 21.0886 74.0000 97.6050 20.1331 76.0000 98.1557 19.1691 78.0000 98.6261 18.1966 80.0000 98.9995 17.2153 82.0000 99.2593 16.2252 84.0000 99.4372 15.2290 86.0000 99.6052 14.2321 88.0000 99.7525 13.2346 90.0000 99.8672 12.2363 92.0000 99.9371 11.2372 94.0000 99.9648 10.2373 96.0000 100.0032 9.2373 98.0000 100.0471 8.2373 100.0000 100.0833 7.2374 102.0000 100.0988 6.2375 104.0000 100.0805 5.2377 106.0000 100.0324 4.2380 108.0000 100.0077 3.2381 110.0000 99.9968 2.2381 112.0000 100.0030
1.2381
北京交通大学海滨学院考试试题 课程名称：数学实验2010-2011第一学期出题教师：数学组适用专业: 09机械, 物流, 土木, 自动化 班级：学号：姓名： 选做题目序号： 1.一对刚出生的幼兔经过一个月可以长成成兔, 成兔再经过一个月后可以 繁殖出一对幼兔. 如果不计算兔子的死亡数, 请用Matlab程序给出在未来24个月中每个月的兔子对数。 解: 由题意每月的成兔与幼兔的数量如下表所示： 1 2 3 4 5 6 ··· 成兔0 1 1 2 3 5··· 幼兔 1 0 1 1 2 3··· 运用Matlab程序： x=zeros(1,24); x(1)=1;x(2)=1； for i=2:24 x(i+1)=x(i)+x(i-1); end x 结果为x = 1 1 2 3 5 8 13 21 3 4 5 5 89 144 233 377 610 987 1597 2584 4181 6765 1094 6 7711 2865 7 46368 2.定积分的过程可以分为分割、求和、取极限三部分, 以1 x e dx 为例, 利用
已学过的Matlab 命令, 通过作图演示计算积分的过程, 并与使用命令int() 直接积分的结果进行比较. 解：根据求积分的过程，我们先对区间[0,1]进行n 等分， 然后针对函数x e 取和，取和的形式为10 1 i n x i e e dx n ξ=≈ ∑ ? ，其中1[ ,]i i i n n ξ-?。这里取i ξ为区间的右端点，则当10n =时，1 x e dx ?可用10 101 1.805610 i i e ==∑ 来近似计算， 当10n =0时，100 100 1 01 =1.7269100 i x i e e dx =≈ ∑?，当10n =000时，10000 10000 1 1 =1.718410000 i x i e e dx =≈ ∑ ?. 示意图如下图，Matlab 命令如下： x=linspace (0,1,21); y=exp(x); y1=y(1:20); s1=sum(y1)/20 y2=y(2:21); s2=sum(y2)/20 plot(x,y); hold on for i=1:20 fill([x(i),x(i+1),x(i+1),x(i),x(i)],[0,0,y(i),y(i),0],'b') end syms k;symsum(exp(k/10)/10,k,1,10)；%n=10 symsum(exp(k/100)/100,k,1,100)；%n=100 symsum(exp(k/10000)/10000,k,1,10000)；%n=10000
大学数学实验 项目一 矩阵运算与方程组求解 实验1 行列式与矩阵 实验目的 掌握矩阵的输入方法. 掌握利用Mathematica (4.0以上版本) 对矩阵进行转置、加、减、数乘、相乘、乘方等运算, 并能求矩阵的逆矩阵和计算方阵的行列式. 基本命令 在Mathematica 中, 向量和矩阵是以表的形式给出的. 1. 表在形式上是用花括号括起来的若干表达式, 表达式之间用逗号隔开. 如输入 {2,4,8,16} {x,x+1,y,Sqrt[2]} 则输入了两个向量. 2. 表的生成函数 (1) 最简单的数值表生成函数Range, 其命令格式如下: Range[正整数n]—生成表{1,2,3,4,…,n }; Range[m, n]—生成表{m ,…,n }; Range[m, n, dx]—生成表{m ,…,n }, 步长为d x . (2) 通用表的生成函数Table. 例如,输入命令 Table[n^3,{n,1,20,2}] 则输出 {1,27,125,343,729,1331,2197,3375,4913,6859} 输入 Table[x*y,{x,3},{y,3}] 则输出 {{1,2,3},{2,4,6},{3,6,9}} 3. 表作为向量和矩阵 一层表在线性代数中表示向量, 二层表表示矩阵. 例如,矩阵 ??? ? ??5432 可以用数表{{2,3},{4,5}}表示. 输入 A={{2,3},{4,5}} 则输出 {{2,3},{4,5}} 命令MatrixForm[A]把矩阵A 显示成通常的矩阵形式. 例如, 输入命令: MatrixForm[A] 则输出 ??? ? ??5432 但要注意, 一般地, MatrixForm[A]代表的矩阵A 不能参与运算. 输入 B={1,3,5,7} 输出为 {1,3,5,7} 输入 MatrixForm[B] 输出为
一牛顿环的各环是否等宽?密度是否均匀?解释原因? 因为环是由空气劈上下表面反射的两束光叠加干涉形成的。劈的上表面变化在横向是不均匀的，故光程差也不是均匀变化的。所以各环是不等宽的环的密度也不是均匀的。各环不等宽，半径小的环宽，越到外边越窄，密度是不均匀的,牛顿环的半径公式是：半径r等于根号下(m+1/2)λR，其中m为环的级数。从公式可以看出，半径和环数并不是线性关系，这样环自然不均匀。计算可以知道，越往外环越密。 二牛顿环的干涉圆环是由哪两束相干光干涉产生的? 半凸透镜下表面和下底面上表面的两束反射光 三电桥由哪几部分组成?电桥平衡的条件? 由电源、开关、检流计桥臂电阻组成。 平衡条件是Rx=(R1/R2)R3 四接通电源后,检流计指针始终向一边偏转,试分析出现这种情况的原因? 指针向一侧偏转就说明发生了电子的定向移动了，这个应该没问题。 指针不偏转，有2种情况吧，其1呢是整个电路发生了断路或其他故障，还1种情况则是流过的电流太小，不足于使电表发生偏转或其偏转的角度肉眼根本看不到。 无论如何调节，检流计指针都不动，电路中可能出现故障是调节臂电阻断路或短路。。无论如何调节，检流计指针始终像一边偏而无法平衡，电路中有可能出现故障是有一个臂（非调节臂）的电阻坏了。（断路或短路） 五什么叫铁磁材料的磁滞现象? 铁磁物质经外磁场磁化到饱和以后，把磁场去掉。这些物质仍保留有剩余磁化强度。需要反方向加磁场才能把这剩余磁化强度变为零。这种现象称为铁磁的磁滞现象。也是说，铁磁材料的磁状态，不仅要看它现在所处的磁场条件；而且还要看它过去的状态。 六如何判断铁磁材料属于软.硬材料? 软磁材料的特点是：磁导率大，矫顽力小，磁滞损耗小，磁滞回线呈长条状；硬磁材料的特点是：剩磁大，矫顽力也大 用光栅方程进行测量的条件是什么? 条件是一束平行光垂直射入光栅平面上，光波发生衍射，即可用光栅方程进行计算。如何实现：使用分光计，光线通过平行光管射入，当狭缝位于透镜的焦平面上时，就能使射在狭缝上的光经过透镜后成为平行光 用光栅方程进行测量,当狭缝太窄或者太宽会怎么样?为什么? 缝太窄，入射光的光强太弱，缝太宽，根据光的空间相干性可以知道，条纹的明暗对比度会下降！ 区别是，太窄了，亮纹会越来越暗，暗纹不变，直到一片黑暗！ 太宽，暗条纹会逐渐加强，明纹不变，直到一片光明！
一、实验内容 P206第六题 function f=wuyan2(c) y=[3.9 5.3 7.2 9.6 12.9 17.1 23.2 31.41 38.6 50.2 62.9 76.0 92.0 106.5 123.2 131.7 150.7 179.3 204.0 226.5 251.4 281.4] t=[0 10 20 30 40 50 60 70 80 90 100 110 120 130 140 150 160 170 180 190 200 210] f=y-c(1)/(1+c(1)/3.9-1)*exp^(-c(2)*t) c0=[1 1] c=lsqnonlin('wuyan2',c0) P206第七题 function f=wuyan1(c) q=[0.4518 0.4862 0.5295 0.5934 0.7171 0.8964 1.0202 1.1963 1.4928 1.6909 1.8548 2.1618 2.6638 3.4634 4.6759 5.8478 6.7885 7.4463 7.8345 8.2068 8.9468 9.7315 10.5172 11.7390 13.6876 ]; k=[0.0911 0.0961 0.1230 0.1430 0.1860 0.2543 0.3121 0.3792 0.4754 0.4410 0.4517 0.5595 0.8080 1.3072 1.7042 2.0019 2.2914 2.4941 2.8406 2.9855 3.2918 3.7214 4.3500 5.5567 7.0477]; l=[4.2361 4.3725 4.5295 4.6436 4.8179 4.9873 5.1282 5.2783 5.4334 5.5329 6.4749 6.5491 6.6152 6.6808 6.7455 6.8065 6.8950 6.9820 7.0637 7.1394 7.2085 7.3025 7.3470 7.4432 7.5200]; f=q-c(1)*k.^c(2).*l.^c(3) c0=[1 1 1] c=lsqnonlin('wuyan1',c0) c = 0.4091 0.6401 1.1446 a=0.4091 α=0.6401 β=1.1446 P239第五题 c=[-20 -30]; A=[1 2;5 4]; b=[20 70]; v1=[0 0]; [x,f,ef,out,lag]=linprog(c,A,b,[],[],v1) z=-f x = 10.0000 5.0000
清华大学数学实验报告4
————————————————————————————————作者: ————————————————————————————————日期： ?
电13 苗键强2011010645
一、实验目的 1.掌握用 MATLAB 软件求解非线性方程和方程组的基本用法, 并对结果作初步分析； 2.练习用非线性方程和方程组建立实际问题的模型并进行求解。 二、实验内容 题目1 【问题描述】 (Q1）小张夫妇以按揭方式贷款买了1套价值２0万元的房子，首付了5万元，每月还款1000元，15年还清。问贷款利率是多少？ (Q2)某人欲贷款50 万元购房,他咨询了两家银行,第一家银行 开出的条件是每月还45０0元,15 年还清；第二家银行开出的条件是每年还４5０00 元,2０年还清。从利率方面看,哪家银行较优惠(简单假设:年利率=月利率×12）? 【分析与解】 假设初始贷款金额为x０,贷款利率为ｐ,每月还款金额为x,第i 个月还完当月贷款后所欠银行的金额为x i,(i=1,2，3,.．....，n）。由题意可知： x1=x0(1+p)?x x2=x0(1+p)2?x(1+p)?x x3=x0(1+p)3?x(1+p)2?x(1+p)?x ……
x n=x0(1+p)n?x(1+p)n?1???x(1+p)?x =x0(1+p)n?x (1+p)n?1 p =0 因而有： x0(1+p)n=x (1+p)n?1 p (1) 则可以根据上述方程描述的函数关系求解相应的变量。 （Ｑ1) 根据公式(１)，可以得到以下方程： 150p(1+p)180?(1+p)180+1=0 设 f(p)=150p(1+p)180?(1+p)180+1，通过计算机程序绘制ｆ（p）的图像以判断解ｐ的大致区间,在Ｍatlab中编程如下： ｆorｉ = 1:25 t = ０．0001＊i; p(i) = t； f(i) =１50*ｔ＊(１+t).^180-(１+t).^180+1; ｅｎd; ｐｌｏｔ(ｐ,ｆ),hｏld ｏn,grid on; 运行以上代码得到如下图像：
第七章 常微分方程初值问题的数值解法 --------学习小结 一、本章学习体会 通过本章的学习，我了解了常微分方程初值问题的计算方法，对于解决那些很难求解出解析表达式的，甚至有解析表达式但是解不出具体的值的常微分方程非常有用。在这一章里求解常微分方程的基本思想是将初值问题进行离散化，然后进行迭代求解。在这里将初值问题离散化的方法有三种，分别是差商代替导数的方法、Taylor 级数法和数值积分法。常微分方程初值问题的数值解法的分类有显示方法和隐式方法，或者可以分为单步法和多步法。在这里单步法是指计算第n+1个y 的值时,只用到前一步的值，而多步法则是指计算第n+1个y 的值时，用到了前几步的值。通过对本章的学习，已经能熟练掌握如何用Taylor 级数法去求解单步法中各方法的公式和截断误差，但是对线性多步法的求解理解不怎么透切，特别是计算过程较复杂的推理。 在本章的学习过程中还遇到不少问题，比如本章知识点多，公式多，在做题时容易混淆，其次对几种R-K 公式的理解不够透彻，处理一个实际问题时，不知道选取哪一种公式，通过课本里面几种方法的计算比较得知其误差并不一样，，这个还需要自己在往后的实际应用中多多实践留意并总结。 二、本章知识梳理 常微分方程初值问题的数值解法一般概念 步长h ，取节点0,(0,1,...,)n t t nh n M =+=，且M t T ≤，则初值问题000 '(,),()y f t y t t T y t y =≤≤?? =?的数值解法的一般形式是 1(,,,...,,)0,(0,1,...,)n n n n k F t y y y h n M k ++==-
常微分方程数值解实验报告 学院:数学与信息科学 专业:信息与计算科学 姓名：郑思义 学号:2０121６52４ 课程：常微分方程数值解 实验一：常微分方程得数值解法 1、分别用Eｕleｒ法、改进得Euler法(预报校正格式)与S—K法求解初值问题。(h＝0、1)并与真解作比较。 1、1实验代码: %欧拉法 functioｎ [ｘ，y]=naeuler(dｙｆｕｎ，xspan,y０,ｈ) %ｄyfun就是常微分方程,xｓpan就是x得取值范围，ｙ０就是初值,h就是步长 x＝ｘsｐan(1）:h:xspan(2)； y(1）=y０； ｆｏrｎ=1:length(x)-１ ｙ（n+１）=ｙ(n)+h*feval（dｙｆun,ｘ(ｎ),ｙ(ｎ））; eｎd ％改进得欧拉法 ｆunction [x,ｍ,ｙ]=naｅuler2(ｄyfun,xspaｎ,y０,h) %ｄyfｕn就是常微分方程,xｓpan就是ｘ得取值范围,y0就是初值，h就是步长。 %返回值ｘ为x取值,ｍ为预报解,ｙ为校正解 x＝xspan(1)：h:xspaｎ(2); y（1)=y0; m=zerｏs（lenｇth(x)－1，１); foｒｎ=1:leｎgｔh(x）-1 k1＝fｅｖal（dｙfun,x(n),y(n）)； y(n＋1)=y（n)+h*k1; ｍ（n)=y（n+1); k2＝ｆevａｌ(dyfuｎ,x(n＋1）,y(n+1）); y（ｎ+1）=y(n)+h*（k１+k2)／2;
end %四阶S—K法 functioｎ [ｘ，ｙ]=rk(ｄyfun,xspａn,y0,h） %dｙｆun就是常微分方程,xｓpaｎ就是ｘ得取值范围,y0就是初值，h就是步长。 ｘ＝xsｐan(1):h：xsｐan(2)； ｙ(１)=y0; for n=1:lｅｎgth(x)-１ ｋ1=fｅｖａl(dyfｕn，x(n），y(ｎ）)； k2=fｅｖal(dyfun,x(ｎ)+ｈ/2，y(ｎ）+(h＊ｋ1)/2); k3=feval(ｄyfuｎ,ｘ(n)+ｈ/２，y(n)+(h*ｋ2)/2); k4=ｆevａl(dｙfun，x(ｎ)+h,y(n)+ｈ＊k3）； y(n+１)＝ｙ(n)+（h/6）*(k1+２*k2+2*k3+k４）; enｄ ％主程序 x=［０:0、1:１］;y＝ｅxp(-ｘ)+x； dyfun=inline('-ｙ+x+1'）； [x1，y1]=naeuler(ｄyｆｕn,[0,１]，1,0、1); ［x2,m，y2]＝naeuler2（dｙfuｎ,[０，1］,1,０、1）; [x3,y3]=rｋ(dyｆuｎ,[0,１］,1,０、１); ploｔ(x，y,'r',x1,y１,'+'，x2,y2,＇*',x3，y３,'o'); xlabｅl('ｘ');yｌabel('y＇); legｅnd(＇y为真解',＇y1为欧拉解＇,'ｙ2为改进欧拉解','y3为S—K解','Locatiｏn','ＮorthWeｓｔ＇)； 1、２实验结果： x 真解y 欧拉解y1 预报值ｍ校正值y2 S—K解y3 0、０１、0000 1、0000 1、０000 １、000０ 0、1 1、0048 1、000０1、００00 １、005０１、0０4８ 0、2 １、0187 １、0100 1、０145 1、０190 1、0187 ０、3 1、0408 １、02９0 1、０３7１1、0412 1、0４0８ ０、4 １、0７03 1、0５６1 1、06７1 1、0708 １、0７03 0、5 1、10６5 1、09０５1、1０3７1、107１1、１0６5 0、6 1、１488 1、１３1４１、1464 １、149４1、１488 0、7 1、1966 １、17８3 1、194５1、19７2 1、１966 0、８1、2493 １、2３05 1、2475 1、2500 1、２４９3 0、9 1、306６1、２87４1、３050 1、3072 １、３066 1、0 1、3679 1、34８７1、36６5 1、３６85 1、36７9
实验一霍尔效应及其应用 【预习思考题】 1．列出计算霍尔系数、载流子浓度n、电导率σ及迁移率μ的计算公式，并注明单位。 霍尔系数，载流子浓度，电导率，迁移率。 2．如已知霍尔样品的工作电流及磁感应强度B的方向，如何判断样品的导电类型？ 以根据右手螺旋定则，从工作电流旋到磁感应强度B确定的方向为正向，若测得的霍尔电压为正，则样品为P型，反之则为N型。 3．本实验为什么要用3个换向开关？ 为了在测量时消除一些霍尔效应的副效应的影响，需要在测量时改变工作电 流及磁感应强度B的方向，因此就需要2个换向开关；除了测量霍尔电压，还要测量A、C间的电位差，这是两个不同的测量位置，又需要1个换向开关。总之，一共需要3个换向开关。 【分析讨论题】 1．若磁感应强度B和霍尔器件平面不完全正交，按式（5.2-5）测出的霍尔系数比实际值大还是小？要准确测定值应怎样进行？ 若磁感应强度B和霍尔器件平面不完全正交，则测出的霍尔系数比实际值偏小。要想准确测定，就需要保证磁感应强度B和霍尔器件平面完全正交，或者设法测量出磁感应强度B和霍尔器件平面的夹角。 2．若已知霍尔器件的性能参数，采用霍尔效应法测量一个未知磁场时，测量误差有哪些来源？ 误差来源有：测量工作电流的电流表的测量误差，测量霍尔器件厚度d的长度测量仪器的测量误差，测量霍尔电压的电压表的测量误差，磁场方向与霍尔器件平面的夹角影响等。 实验二声速的测量 【预习思考题】 1. 如何调节和判断测量系统是否处于共振状态？为什么要在系统处于共振的条件下进行声速测定？ 答：缓慢调节声速测试仪信号源面板上的“信号频率”旋钮，使交流毫伏表指针指示达到最大（或晶体管电压表的示值达到最大），此时系统处于共振状态，显示共振发生的信号指示灯亮，信号源面板上频率显示窗口显示共振频率。在进行声速测定时需要测定驻波波节的位置，当发射换能器S1处于共振状态时，发射的超声波能量最大。若在这样一个最佳状态移动S1至每一个波节处，媒质压缩形变最大，则产生的声压最大，接收换能器S2接收到的声压为最大，转变成电信号，晶体管电压表会显示出最大值。由数显表头读出每一个电压最大值时的位置，即对应的波节位置。因此在系统处于共振的条件下进行声速测定，可以容易和准确地测定波节的位置，提高测量的准确度。 2. 压电陶瓷超声换能器是怎样实现机械信号和电信号之间的相互转换的？ 答：压电陶瓷超声换能器的重要组成部分是压电陶瓷环。压电陶瓷环由多晶结构的压电材料制成。这种材料在受到机械应力，发生机械形变时，会发生极化，同时在极化方向产生电场，这种特性称为压电效应。反之，如果在压电材料上加交
实验2 方程模型及其求解算法 一、实验目的及意义 [1] 复习求解方程及方程组的基本原理和方法； [2] 掌握迭代算法； [3] 熟悉MATLAB软件编程环境；掌握MATLAB编程语句(特别是循环、条件、控制等语句)； [4] 通过范例展现求解实际问题的初步建模过程； 通过该实验的学习，复习和归纳方程求解或方程组求解的各种数值解法（简单迭代法、二分法、牛顿法、割线法等），初步了解数学建模过程。这对于学生深入理解数学概念，掌握数学的思维方法，熟悉处理大量的工程计算问题的方法具有十分重要的意义。 二、实验内容 1．方程求解和方程组的各种数值解法练习 2．直接使用MATLAB命令对方程和方程组进行求解练习 3．针对实际问题，试建立数学模型，并求解。 三、实验步骤 1．开启软件平台——MATLAB，开启MATLAB编辑窗口； 2．根据各种数值解法步骤编写M文件 3．保存文件并运行； 4．观察运行结果(数值或图形)； 5．根据观察到的结果写出实验报告，并浅谈学习心得体会。 四、实验要求与任务 基础实验 1．用图形放大法求解方程x sin(x) = 1. 并观察该方程有多少个根。 画出图形程序： x=-10:0.01:10; y=x.*sin(x)-1; y1=zeros(size(x)); plot(x,y,x,y1) MATLAB运行结果：
-10-8-6-4-20246810 -8-6 -4 -2 2 4 6 8 扩大区间画图程序： x=-50:0.01:50; y=x.*sin(x)-1; y1=zeros(size(x)); plot(x,y,x,y1) MATLAB 运行结果： -50-40-30-20-1001020304050 由上图可知，该方程有偶数个无数的根。
习 题(参考答案) 2．指出下列测量值为几位有效数字，哪些数字是可疑数字，并计算相对不确定度。 (1) g ＝(9.794±0.003)m ·s 2 - 答：四位有效数字，最后一位“4”是可疑数字，%031.0%100794 .9003 .0≈?= gr U ； (2) e ＝(1.61210±0.00007)?10 19 - C 答：六位有效数字，最后一位“0”是可疑数字，%0043.0%10061210 .100007 .0≈?= er U ； (3) m ＝(9.10091±0.00004) ?10 31 -kg 答：六位有效数字，最后一位“1”是可疑数字，%00044.0%10010091 .900004 .0≈?= mr U ； (4) C ＝(2.9979245±0.0000003)8 10?m/s 答：八位有效数字，最后一位“5”是可疑数字 1．仪器误差为0.005mm 的螺旋测微计测量一根直径为D 的钢丝，直径的10次测量值如下表： 试计算直径的平均值、不确定度(用D 表示)和相对不确定度(用Dr 表示)，并用标准形式表示测量结果。 解： 平均值 mm D D i i 054.210110 1 ==∑=
标准偏差： mm D D i i D 0029.01 10)(10 1 2 ≈--= ∑=σ 算术平均误差： m m D D i i D 0024.010 10 1 ≈-= ∑=δ 不确定度A 类分量mm U D A 0029.0==σ， 不确定度B 类分量mm U B 005.0=?=仪 ∴ 不确定度mm U U U B A D 006.0005.00029.0222 2≈+=+= 相对不确定度%29.0%100054 .2006 .0%100≈?=?= D U U D Dr 钢丝的直径为：%29.0)006.0054.2(=±=Dr D mm D 或 不确定度A 类分量mm U D A 0024.0==δ ， 不确定度B 类分量mm U B 005.0=?=仪 ∴ 不确定度mm U U U B A D 006.0005.00024.0222 2≈+=+= 相对不确定度%29.0%100054 .2006 .0%100≈?=?= D U U D Dr 钢丝的直径为： %29.0)006.0054.2(=±=Dr D mm D ，%00001.0%1009979245 .20000003 .0≈?= Cr U 。 3．正确写出下列表达式 （1）km km L 310)1.01.3()1003073(?±=±= （2）kg kg M 4 10)01.064.5()13056430(?±=±= （3）kg kg M 4 10)03.032.6()0000030.00006320.0(-?±=±= （4）s m s m V /)008.0874.9(/)00834 .0873657.9(±=±= 4．试求下列间接测量值的不确定度和相对不确定度，并把答案写成标准形式。
东华大学M A T L A B数学实验第二版答案(胡良 剑) -CAL-FENGHAI-(2020YEAR-YICAI)_JINGBIAN
数学实验答案 Chapter 1 Page20,ex1 (5) 等于[exp(1),exp(2);exp(3),exp(4)] (7) 3=1*3, 8=2*4 (8) a为各列最小值，b为最小值所在的行号 (10) 1>=4,false, 2>=3,false, 3>=2, ture, 4>=1,ture (11) 答案表明：编址第2元素满足不等式(30>=20)和编址第4元素满足不等式(40>=10) (12) 答案表明：编址第2行第1列元素满足不等式(30>=20)和编址第2行第2列元素满足不等式(40>=10) Page20, ex2 (1)a, b, c的值尽管都是1，但数据类型分别为数值，字符，逻辑，注意a与c 相等，但他们不等于b (2)double(fun)输出的分别是字符a,b,s,(,x,)的ASCII码 Page20,ex3 >> r=2;p=0.5;n=12; >> T=log(r)/n/log(1+0.01*p) Page20,ex4 >> x=-2:0.05:2;f=x.^4-2.^x; >> [fmin,min_index]=min(f) 最小值最小值点编址 >> x(min_index) ans = 0.6500 最小值点 >> [f1,x1_index]=min(abs(f)) 求近似根--绝对值最小的点 f1 = 0.0328 x1_index = 24 >> x(x1_index) ans = -0.8500 >> x(x1_index)=[];f=x.^4-2.^x; 删去绝对值最小的点以求函数绝对值次小的点>> [f2,x2_index]=min(abs(f)) 求另一近似根--函数绝对值次小的点 f2 = 0.0630 x2_index = 65 >> x(x2_index) ans =
大学数学实验心得体会 [模版仅供参考，切勿通篇使用] 大学数学实验心得体会（一） 数学，在整个人类生命进程中至关重要，从小学到中学，再到大学，乃至更高层次的科学研究都离不开数学，随着时代的发展，人们越来越重视数学知识的应用，对数学课程提出了更高层次的要求，于是便诞生了数学实验。 学期最初，大学数学实验对于我们来说既熟悉又陌生，在我们的记忆中，我们做过物理实验、化学实验、生物实验，故然我们以为数学实验与它们一样，当我们在网上搜索有关数学实验的信息时，我们才知道，大学数学实验作为一门新兴的数学课程在近十年来取得了迅速的发展。数学实验以计算机技术和数学软件为载体，将数学建模的思想和方法融入其中，现在已经成为一种潮流。 当我们怀着好奇的心情走进屈静国老师的数学实验课堂时，我们才渐渐懂得，数学实验是一门有关计算机软件的课程，就像c语言一样，需要编辑运行程序，从而进行数学运算，它不需要自己来运算，就像计算器一样，只要我们自己记下重要程序语句，输入运行程序，便可得到运行结果，大大降低了我们的运算量，
给我们生活带来许多便捷，在大一时，我学过c语言，由于这样的基础，让我能够更快的学会并应用此软件。 时间飞逝，转眼间，我们就要结课了，这学期我们学习了mathematics的基础，微积分实验，线性代数实验，概率论与数理统计实验，数值计算方法及实验。通过这学期的学习，我也积累了些自己的学习方法和心得。首先，我们要在平时上课牢记那些mathematics语言和公式，那些东西就想单词和公式一样，只需要背诵;然后，我们要看几遍书，并多看一下例题;最后，我们要多应用mathematics软件去练习。正所谓熟能生巧，我坚信，只要我们能够做到这三步，我们就能很好的掌握这门课程。 通过学习使用数学软件，数学实验建模，使我们能够从实际问题出发，认真分析研究，建立简单数学模型，然后借助先进的计算机技术，最终找出解决实际问题的一种或多种方案，从而提高了我们的数学思维能力，为我们参加数学竞赛和数学建模打下了坚实的基础，同时也为我们进一步深造和参加工作打下一定的实践基础! 大学数学实验心得体会（二） 在此期间我充分利用研修活动时间学习，感到既有辛苦，又有收获。既有付出，又有新所得。这次远程研修让我有幸与专家和各地的数学精英们交流，面对每次探讨的主题，大家畅所欲言，
第四章常微分方程数值解 [课时安排]6学时 [教学课型]理论课 [教学目的和要求] 了解常微分方程初值问题数值解法的一些基本概念，如单步法和多步法，显式和隐式，方法的阶数，整体截断误差和局部截断误差的区别和关系等；掌握一阶常微分方程初值问题的一些常用的数值计算方法，例如欧拉（Euler）方法、改进的欧拉方法、龙贝－库塔（Runge-Kutta）方法、阿达姆斯（Adams）方法等，要注意各方法的特点及有关的理论分析；掌握构造常微分方程数值解的数值积分的构造方法和泰勒展开的构造方法的基本思想，并能具体应用它们导出一些常用的数值计算公式及评估截断误差；熟练掌握龙格－库塔（Ｒ－Ｋ）方法的基本思想，公式的推导，Ｒ－Ｋ公式中系数的确定，特别是能应用“标准四阶Ｒ－Ｋ公式”解题；掌握数值方法的收敛性和稳定性的概念，并能确定给定方法的绝对稳定性区域。[教学重点与难点] 重点：欧拉方法，改进的欧拉方法，龙贝－库塔方法。 难点：R—K方法，预估-校正公式。 [教学内容与过程] 4.1 引言 本章讨论常微分方程初值问题 (4.1.1) 的数值解法，这也是科学与工程计算经常遇到的问题，由于只有很特殊的方程能用解析方法求解，而用计算机求解常微分方程的初值问题都要采用数值方法.通常我们假定(4.1.1)中 f(x,y)对y满足Lipschitz条件，即存在常数L＞0，使对，有 (4.1.2) 则初值问题(4.1.1)的解存在唯一. 假定(4.1.1)的精确解为，求它的数值解就是要在区间上的一组离散点 上求的近似.通常取 ，h称为步长，求(4.1.1)的数值解是按节点的顺序逐步 推进求得.首先，要对方程做离散逼近，求出数值解的公式，再研究公式的局部截
1. 分别用Euler 法和ode45解下列常微分方程并与解析解比较： (1) 30,1)0(,<<=+='x y y x y 编写Euler 法的matlab 函数，命名为euler.m function [t,y]=euler(odefun,tspan,y0,h)t=tspan(1):h:tspan(2);y(1)=y0; for i=1:length(t)-1 y(i+1)=y(i)+h*feval(odefun,t(i),y(i));end t=t';y=y'; 下面比较三者的差别：% ode45 odefun=inline('x+y','x','y');[x1,y1]=ode45(odefun,[0,3],1);plot(x1,y1,'ko');pause hold on ;% Euler·¨ [x2,y2]=euler(odefun,[0,3],1,0.05);plot(x2,y2,'r+');pause hold on ; % 解析解 y0=dsolve('Dy=t+y','y(0)=1');ezplot(y0,[0,3]);pause hold off ; legend('ode45','euler 法','解析解');
Euler 法只有一阶精度，所以实际应用效率比较差，而ode45的效果比较好，很接近真实值。 (2) 2 0.01()2sin(),(0)0,(0)1,05 y y y t y y t ''''-+===<<先写M 文件ex1_2fun.m function f=ex1_2fun(t,y)f(1)=y(2); f(2)=0.01*y(2).^2-2*y(1)+sin(t);f=f(:);% ode45 [t1,y1]=ode45(@ex1_2fun,[0,5],[0;1]);plot(t1,y1(:,1),'ko'); % 解析解 s=dsolve('D2y-0.01*(Dy)^2+2*y=sin(t)','y(0)=0','Dy(0)=1','t') s = [ empty sym ] %由此可知该微分方程无解析解2. 求一通过原点的曲线，它在处的切线斜率等于若上限增为1.58，1.60会(,)x y 2 2,0 1.57.x y x +< 大学物理实验课后答 案 （1）利用f=(D+d)(D-d)/4D 测量凸透镜焦距有什么优点？ 答这种方法可以避免透镜光心位置的不确定而带来的测量物距和像距的误差。 （2）为什么在本实验中利用1/u+1/v=1/f 测焦距时，测量u和v都用毫米刻度的米尺就可以满足要求？设透镜由于色差和非近轴光线引起的误差是 1％。 答设物距为20cm，毫米刻度尺带来的最大误差为0.5mm，其相对误差为 0.25％，故没必要用更高精度的仪器。 （3）如果测得多组u，v值，然后以u+v为纵轴，以uv为横轴，作出实验的曲线属于什么类型，如何利用曲线求出透镜的焦距f。 答直线；1/f为直线的斜率。 （4）试证：在位移法中，为什么物屏与像屏的间距D要略大于4f？ 由f=(D+d)(D-d)/4D → D2-4Df=d2→ D(D-4f)=d2 因为d>0 and D>0 故D>4f 1.避免测量u、ν的值时，难于找准透镜光心位置所造成的误差。 2.因为实验中，侧的值u、ν、f都相对较大，为十几厘米到几十厘米左右，而误差为1%，即一毫米到几毫米之间，所以可以满足要求。 3.曲线为曲线型曲线。透镜的焦距为基斜率的倒数。 ①当缝宽增加一倍时,衍射光样的光强和条纹宽度将会怎样变化?如缝宽减半,又怎样改变? 答: a增大一倍时, 光强度↑；由a=Lλ/b ，b减小一半 a减小一半时, 光强度↓；由a=Lλ/b ，b增大一倍。 ②激光输出的光强如有变动,对单缝衍射图象和光强分布曲线有无影响?有何影响? 答:由b=Lλ/a.无论光强如何变化,只要缝宽不变,L不变,则衍射图象的光强分布曲线不变 (条纹间距b不变)；整体光强度↑或者↓。 ③用实验中所应用的方法是否可测量细丝直径？其原理和方法如何？ 答：可以，原理和方法与测单狭缝同。 ④本实验中，λ=632。8nm,缝宽约为5*10^-3㎝,屏距L为50㎝。试验证： 是否满足夫朗和费衍射条件？ 答：依题意： Lλ=（50*10^-2）*（632.8*10^-9）=3.164*10^-7 a^2/8=(5*10^-5)^2/8=3.1*10^-10 所以Lλ< 第七章 常微分方程数值解法 常微分方程中只有一些典型方程能求出初等解(用初等函数表示的解)，大部分的方程是求不出初等解的。另外，有些初值问题虽然有初等解，但由于形式太复杂不便于应用。因此，有必要探讨常微分方程初值问题的数值解法。本章主要介绍一阶常微分方程初值问题的欧拉法、龙格-库塔法、阿达姆斯方法，在此基础上推出一阶微分方程组与高阶方程初值问题的 数值解法；此外，还将简要介绍求解二阶常微分方程值问题的差分方法、试射法。 第一节 欧拉法 求解常微分方程初值问题 ?????==0 0)() ,(y x y y x f dx dy （1） 的数值解，就是寻求准确解)(x y 在一系列离散节点 <<<< 常微分方程的边值问题和本征值问题 一、问题描述 利用搜索法和弦割法，得到该常微分方程的本征值，再利用打靶法计算多个本征值。 二、解决方法 （一）搜索法 1.先随便猜测k的一个试验值,程序中令k=1 2.由Numerov算法 根据本题的条件，kn+1=kn=kn-1=k,s=0,得到yn+2,yn+1,yn间的迭代公式 令con=(k*h)^2/12 yn+2=2*(1-5*con)*yn+1/(1+con)-yn 3自己给定φ的初始条件，然后利用公式得到边界值φ(1) 4.然后以小的步长dk增加k值，这里令dk=1,每当φ(1)改变符号时，就将步长减半后倒退回来重复 5.当步长小于所要求的容许误差时终止程序，此时的k值即为所求。 （二）弦割法 1.随便猜测两个k值，这里令k0=1,k1=2 2.自己给定φ的初始条件，对两个k值分别利用上述公式进行迭代，得到边界值y1(1)和y2(1)。 3.比较y1(1)和y2(1)的绝对值大小。若绝对值大，说明对应的k值距离本征值距离较远。 4.将（k0+k1）/2赋给k2,边界值绝对值小的对应的k值保持不变，边界值绝对值大的对应k值重新定位k2的值。 5.重复进行实验，当y1(1)和y(2)的差的绝对值小于容许误差时终止程序。此时k1的值即为所求。 当搜索法和弦割法大致求出了一个本征值后，利用打靶法，调整k值再度进行搜索，得到多个本征值，绘出其中一个本征值对应的函数图像，观察其性质。 三、程序实现 1.搜索法 subroutine add(t,y0,y1) !利用子程序表示函数值的迭代 implicit none real(8)::t,h,con,y0,y1,y2 （1）利用f=(D+d)(D-d)/4D 测量凸透镜焦距有什么优点？ 答这种方法可以避免透镜光心位置的不确定而带来的测量物距和像距的误差。 （2）为什么在本实验中利用1/u+1/v=1/f 测焦距时，测量u和v都用毫米刻度的米尺就可以满足要求？设透镜由于色差和非近轴光线引起的误差是1％。 答设物距为20cm，毫米刻度尺带来的最大误差为0.5mm，其相对误差为 0.25％，故没必要用更高精度的仪器。 （3）如果测得多组u，v值，然后以u+v为纵轴，以uv为横轴，作出实验的曲线属于什么类型，如何利用曲线求出透镜的焦距f。 答直线；1/f为直线的斜率。 （4）试证：在位移法中，为什么物屏与像屏的间距D要略大于4f？ 由f=(D+d)(D-d)/4D → D2-4Df=d2→ D(D-4f)=d2 因为d>0 and D>0 故D>4f 1.避免测量u、ν的值时，难于找准透镜光心位置所造成的误差。 2.因为实验中，侧的值u、ν、f都相对较大，为十几厘米到几十厘米左右，而误差为1%，即一毫米到几毫米之间，所以可以满足要求。 3.曲线为曲线型曲线。透镜的焦距为基斜率的倒数。 ①当缝宽增加一倍时,衍射光样的光强和条纹宽度将会怎样变化?如缝宽减半,又怎样改变? 答: a增大一倍时, 光强度↑；由a=Lλ/b ，b减小一半 a减小一半时, 光强度↓；由a=Lλ/b ，b增大一倍。 ②激光输出的光强如有变动,对单缝衍射图象和光强分布曲线有无影响?有何影响? 答:由b=Lλ/a.无论光强如何变化,只要缝宽不变,L不变,则衍射图象的光强分布曲线不变 (条纹间距b不变)；整体光强度↑或者↓。 ③用实验中所应用的方法是否可测量细丝直径？其原理和方法如何？ 答：可以，原理和方法与测单狭缝同。 ④本实验中，λ=632。8nm,缝宽约为5*10^-3㎝,屏距L为50㎝。试验证： 是否满足夫朗和费衍射条件？ 答：依题意： Lλ=（50*10^-2）*（632.8*10^-9）=3.164*10^-7 a^2/8=(5*10^-5)^2/8=3.1*10^-10 所以Lλ< 二、判断题（“对”在题号前（）中打√×）（10分） （√）1、误差是指测量值与真值之差，即误差=测量值－真值，如此定义的误差反映的是测量值偏离真值的大小和方向，既有大小又有正负符号。 （×）2、残差（偏差）是指测量值与其算术平均值之差，它与误差定义一样。（√）3、精密度是指重复测量所得结果相互接近程度，反映的是随机误差大小的程度。 （√）4、测量不确定度是评价测量质量的一个重要指标，是指测量误差可能出现的范围。 （×）7、分光计设计了两个角游标是为了消除视差。 （×）9、调节气垫导轨水平时发现在滑块运动方向上不水平，应该先调节单脚螺钉再调节双脚螺钉。 （×）10、用一级千分尺测量某一长度（Δ仪=0.004mm），单次测量结果为N=8.000mm，用不确定度评定测量结果为N=（8.000±0.004）mm。 三、简答题（共15分） 1．示波器实验中，（1）CH1（x）输入信号频率为50Hz，CH2（y）输入信号频率为100Hz；（2）CH1（x）输入信号频率为150Hz，CH2（y）输入信号频率为50Hz；画出这两种情况下，示波器上显示的李萨如图形。（8分） 差法处理数据的优点是什么？（7分） 答：自变量应满足等间距变化的要求，且满足分组要求。（4分） 优点：充分利用数据；消除部分定值系统误差 四、计算题（20分，每题10分） 1、用1/50游标卡尺，测得某金属板的长和宽数据如下表所示，求金属板的面 解：（1）金属块长度平均值：)(02.10mm L = 长度不确定度： )(01.03/02.0mm u L == 金属块长度为：mm L 01.002.10±= %10.0=B （2分） （2）金属块宽度平均值：)(05.4mm d = 宽度不确定度： )(01.03/02.0mm u d == 金属块宽度是：mm d 01.005.4±= %20.0=B （2分） （3）面积最佳估计值：258.40mm d L S =?= 不确定度：2222222 221.0mm L d d s L s d L d L S =+=??? ????+??? ????=σσσσσ 相对百分误差：B ＝%100?S s σ=0.25％ （4分） （4）结果表达：21.06.40mm S ±= B ＝0.25％ （2分） 注：注意有效数字位数，有误者酌情扣 5、测量中的千分尺的零点误差属于已定系统误差；米尺刻度不均匀的误差属于未大学物理实验课后答案教学内容
常微分方程数值解法
(完整版)计算物理实验报告常微分方程
大学物理实验课后答案
《大学物理实验》模拟试卷与答案
相关主题
文本预览