数值分析经典例题
1. y' = y , x [0,1] ,y (0) =1 , h = 0.1。
1求解析解。
2 Eular法
3 R-K法
1解析法 ○
在MATLAB 命令窗口执行
clear
>> x=0:0.1:1;
>> y=exp(x);
>> c=[y]'
c =
1.[**************]
1.[**************]
1.[**************]
1.[**************]
1.[**************]
1.[**************]
1.[**************]
2.[**************]
2.[**************]
2.[**************]
2.[**************]
2Euler 法 ○
在Matlab 中建立M 文件如下：
function [x,y]=euler1(dyfun,xspan,y0,h)
x=xspan(1):h:xspan(2);y(1)=y0;
for n=1:length(x)-1
y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));
end
x=x';y=y'
在MATLAB 命令窗口执行
clear
>> dyfun=inline('y+0*x');
>> [x,y]=euler1(dyfun,[0,1],1,0.1);
>> [x,y]
得到
ans =
0 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 2.[**************]
0.[**************] 2.[**************]
1.[**************] 2.[**************]
3R-K 法（龙格-库塔法） ○
在本题求解中，采用经典4阶龙格-库塔法
首先在Matlab 的M 文件窗口对4阶龙格-库塔算法进行编程：
function [x,y]=RungKutta41(dyfun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
k1=h*feval(dyfun,x(n),y(n));
k2=h*feval(dyfun,x(n)+h/2,y(n)+1/2*k1);
k3=h*feval(dyfun,x(n)+h/2,y(n)+1/2*k2);
k4=h*feval(dyfun,x(n+1)+h,y(n)+k3);
y(n+1)=y(n)+(k1+2*k2+2*k3+k4)/6;
end
在MATLAB 命令窗口执行
clear
>> dyfun=inline('y','x','y');
>> [x,y]=RungKutta41(dyfun,0,1,0.1,10);
>> c=[x;y]'
得到
c =
0 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 2.[**************]
0.[**************] 2.[**************]
0.[**************] 2.[**************]
1.[**************] 2.[**************]
4绘图 ○
' 解析法','Euler 法','R-K 法' 绘制如下
在MATLAB 命令窗口执行
clear
>> x=0:0.1:1;
>> y1=exp(x);
>> dyfun=inline('y+0*x');
>> [x,y2]=euler1(dyfun,[0,1],1,0.1);
>> dyfun=inline('y','x','y');
>> [x,y3]=RungKutta41(dyfun,0,1,0.1,10);
>> plot(x,y1,'*')
hold on
plot(x,y2,'g','LineWidth',2)
plot(x,y3,'b','LineWidth',2)
legend('解析法','Euler 法','R-K 法')
2. 一个具有1400kg 初始重量的小火箭，带有1040kg 的燃料，点燃后垂直向上运动，火箭内的燃料以18kg/s的速率燃烧，提供31000N 的推力。试求燃料燃尽之前火箭的运动规律。（高度h ，速度v ，加速度a 随时间变化的情况）假定空气阻
22力D ∝ v , D=k*v, k=0.3822 N*S2/m2 。
解：
本题利用Matlab 求解过程如下：
在MATLAB 命令窗口执行
>> clear
>> a=0;
>> v=0;
>> h=0;
>> aa(1:59)=0;
>> vv(1:59)=0;
>> hh(1:59)=0;
>> for t=0:58
aa(t+1)=a;
vv(t+1)=v;
hh(t+1)=h;
l1=(9.8*18*t-0.3822*v^2+31000-1400*9.8)/( 1400-18*t);
l2=(9.8*18*(t+0.5)-0.3822*(v+0.5*l1)^2+31000-1400*9.8)/(1400-18*(t+0.5)); l3=(9.8*18*(t+0.5)-0.3822*(v+0.5*l2)^2+31000-1400*9.8)/(1400-18*(t+0.5)); l4=(9.8*18*(t+1)-0.3822*(v+l3)^2+31000-1400*9.8)/( 1400-18*(t+1)); h=h+v+(l1+l2+l3)/6;
v=v+(l1+2*l2+2*l3+l4)/6;
a=l1;
end
>> t=0:1:58;
○1输入命令
>> plot(t, aa,'r','LineWidth',2)
得燃料燃尽前火箭飞行加速度变化曲线
运行
>> [t;aa]'
得火箭飞行中加速度具体数据
ans =
0 0
1.[**************] 12.[**************]
2.[**************] 12.[**************]
3.[**************] 12.[**************]
4.[**************] 12.[**************]
5.[**************] 12.[**************]
6.[**************] 12.[**************]
7.[**************] 12.[**************]
8.[**************] 12.[**************]
9.[**************] 11.[**************]
10.[**************] 11.[**************]
11.[**************] 10.[**************]
12.[**************] 10.[**************]
13.[**************] 9.[**************]
14.[**************] 9.[**************]
15.[**************] 8.[**************]
16.[**************] 7.[**************]
17.[**************]
7.[**************]
18.[**************] 6.[**************]
19.[**************] 6.[**************]
20.[**************] 5.[**************]
21.[**************] 5.[**************]
22.[**************] 4.[**************]
23.[**************] 4.[**************]
24.[**************] 3.[**************]
25.[**************] 3.[**************]
26.[**************] 3.[**************]
27.[**************]
28.[**************]
29.[**************]
30.[**************]
31.[**************]
32.[**************]
33.[**************]
34.[**************]
35.[**************]
36.[**************]
37.[**************]
38.[**************]
39.[**************]
40.[**************]
41.[**************]
42.[**************]
43.[**************]
44.[**************]
45.[**************]
46.[**************]
47.[**************]
48.[**************]
49.[**************]
50.[**************]
51.[**************]
52.[**************]
53.[**************]
54.[**************]
55.[**************]
56.[**************]
57.[**************]
58.[**************]
○2输入命令 2.[**************] 2.[**************] 2.[**************] 2.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************]
plot(t, vv,'g','LineWidth',2)
得燃料燃尽前火箭飞行速度变化曲线
运行
>> [t;vv]'
得火箭飞行中速度具体数据
ans =
0 0 1 12.[1**********]38 2 25.[1**********]97 3 37.[1**********]46 4 50.[1**********]84 5 63.[1**********]58 6 76.[1**********]68 7 88.[1**********]62 8 100.[1**********] 9 112.[1**********]8 10 123.[1**********]1 11 133.[1**********] 12
143.[1**********]3
13 153.[1**********]5 14 161.[1**********]2 15 170.[1**********]3 16 177.[1**********]6 17 184.[1**********]9 18 190.[1**********]2 19 196.[1**********]1 20 201.[1**********]5 21 206.[1**********]1 22 211.[1**********]1 23 214.[1**********]3 24 218.[1**********]7 25 221.[1**********]9 26 224.[1**********]6 27 227.[1**********]2 28 229.[1**********]4 29 231.[1**********]5 30 233.[1**********]9 31 235.[1**********]9 32 237.[1**********]1 33 239.[1**********]9 34 240.[1**********]7 35 242.[1**********]8 36 243.[1**********]4 37 244.[1**********]1 38 245.[1**********]5 39 247.[1**********]2 40 248.[1**********]4 41 249.[1**********]7 42 250.[1**********]8 43 251.[1**********]2 44 252.[1**********]3 45 253.[1**********]3 46 254.[1**********]2 47 255.[1**********]5 48 256.[1**********]4 49 257.[1**********]7 50 258.[1**********]6 51 259.[1**********]1 52 260.[1**********]9 53 261.[1**********] 54 262.[1**********]2 55 263.[1**********]7 56 264.[1**********]6
57 265.[1**********]2 58 266.[1**********]1 3输入命令 ○
>> plot(t, hh,'b','LineWidth',2)
得燃料燃尽前火箭飞行高度时间变化曲线
运行
>> [t;hh]'
得火箭飞行过程中高度具体数据
ans =
0 0 1 6.[1**********]771 2 25.[1**********]62 3 56.[1**********]49 4 100.[1**********]5 5 158.[1**********]4 6 227.[1**********]1 7 310.[1**********] 8
404.[1**********]7
10 628.[1**********]6 11 757.[1**********]4 12 895.[1**********]7 13 1044.[1**********] 14 1202.[1**********] 15 1368.[1**********] 16 1542.[1**********] 17 1723.1584828712 18 1910.[1**********] 19 2104.[1**********] 20 2304.[1**********] 21 2508.[1**********] 22 2717.[1**********] 23 2930.[1**********] 24 3147.[1**********] 25 3367.[1**********] 26 3590.[1**********] 27 3816.[1**********] 28 4045.[1**********] 29 4276.[1**********] 30 4509.[1**********] 31 4743.[1**********] 32 4980.[1**********] 33 5218.[1**********] 34 5458.[1**********] 35 5699.[1**********] 36 5942.[1**********] 37 6186.[1**********] 38 6431.[1**********] 39 6678.[1**********] 40 6926.[1**********] 41 7174.[1**********] 42 7424.[1**********] 43 7675.[1**********] 44 7928.[1**********] 45 8181.[1**********] 46 8435.508981705 47 8690.[1**********] 48 8947.[1**********] 49 9204.[1**********] 50 9462.[1**********] 51 9721.[1**********] 52 9982.2660536508
54 10505.8206883216 55 10769.0654797243 56 11033.2819028965 57 11298.4652041564 58 11564.6107884638
数值分析经典例题
1. y' = y , x [0,1] ,y (0) =1 , h = 0.1。
1求解析解。
2 Eular法
3 R-K法
1解析法 ○
在MATLAB 命令窗口执行
clear
>> x=0:0.1:1;
>> y=exp(x);
>> c=[y]'
c =
1.[**************]
1.[**************]
1.[**************]
1.[**************]
1.[**************]
1.[**************]
1.[**************]
2.[**************]
2.[**************]
2.[**************]
2.[**************]
2Euler 法 ○
在Matlab 中建立M 文件如下：
function [x,y]=euler1(dyfun,xspan,y0,h)
x=xspan(1):h:xspan(2);y(1)=y0;
for n=1:length(x)-1
y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));
end
x=x';y=y'
在MATLAB 命令窗口执行
clear
>> dyfun=inline('y+0*x');
>> [x,y]=euler1(dyfun,[0,1],1,0.1);
>> [x,y]
得到
ans =
0 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 2.[**************]
0.[**************] 2.[**************]
1.[**************] 2.[**************]
3R-K 法（龙格-库塔法） ○
在本题求解中，采用经典4阶龙格-库塔法
首先在Matlab 的M 文件窗口对4阶龙格-库塔算法进行编程：
function [x,y]=RungKutta41(dyfun,x0,y0,h,N)
x=zeros(1,N+1);y=zeros(1,N+1);x(1)=x0;y(1)=y0;
for n=1:N
x(n+1)=x(n)+h;
k1=h*feval(dyfun,x(n),y(n));
k2=h*feval(dyfun,x(n)+h/2,y(n)+1/2*k1);
k3=h*feval(dyfun,x(n)+h/2,y(n)+1/2*k2);
k4=h*feval(dyfun,x(n+1)+h,y(n)+k3);
y(n+1)=y(n)+(k1+2*k2+2*k3+k4)/6;
end
在MATLAB 命令窗口执行
clear
>> dyfun=inline('y','x','y');
>> [x,y]=RungKutta41(dyfun,0,1,0.1,10);
>> c=[x;y]'
得到
c =
0 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 1.[**************]
0.[**************] 2.[**************]
0.[**************] 2.[**************]
0.[**************] 2.[**************]
1.[**************] 2.[**************]
4绘图 ○
' 解析法','Euler 法','R-K 法' 绘制如下
在MATLAB 命令窗口执行
clear
>> x=0:0.1:1;
>> y1=exp(x);
>> dyfun=inline('y+0*x');
>> [x,y2]=euler1(dyfun,[0,1],1,0.1);
>> dyfun=inline('y','x','y');
>> [x,y3]=RungKutta41(dyfun,0,1,0.1,10);
>> plot(x,y1,'*')
hold on
plot(x,y2,'g','LineWidth',2)
plot(x,y3,'b','LineWidth',2)
legend('解析法','Euler 法','R-K 法')
2. 一个具有1400kg 初始重量的小火箭，带有1040kg 的燃料，点燃后垂直向上运动，火箭内的燃料以18kg/s的速率燃烧，提供31000N 的推力。试求燃料燃尽之前火箭的运动规律。（高度h ，速度v ，加速度a 随时间变化的情况）假定空气阻
22力D ∝ v , D=k*v, k=0.3822 N*S2/m2 。
解：
本题利用Matlab 求解过程如下：
在MATLAB 命令窗口执行
>> clear
>> a=0;
>> v=0;
>> h=0;
>> aa(1:59)=0;
>> vv(1:59)=0;
>> hh(1:59)=0;
>> for t=0:58
aa(t+1)=a;
vv(t+1)=v;
hh(t+1)=h;
l1=(9.8*18*t-0.3822*v^2+31000-1400*9.8)/( 1400-18*t);
l2=(9.8*18*(t+0.5)-0.3822*(v+0.5*l1)^2+31000-1400*9.8)/(1400-18*(t+0.5)); l3=(9.8*18*(t+0.5)-0.3822*(v+0.5*l2)^2+31000-1400*9.8)/(1400-18*(t+0.5)); l4=(9.8*18*(t+1)-0.3822*(v+l3)^2+31000-1400*9.8)/( 1400-18*(t+1)); h=h+v+(l1+l2+l3)/6;
v=v+(l1+2*l2+2*l3+l4)/6;
a=l1;
end
>> t=0:1:58;
○1输入命令
>> plot(t, aa,'r','LineWidth',2)
得燃料燃尽前火箭飞行加速度变化曲线
运行
>> [t;aa]'
得火箭飞行中加速度具体数据
ans =
0 0
1.[**************] 12.[**************]
2.[**************] 12.[**************]
3.[**************] 12.[**************]
4.[**************] 12.[**************]
5.[**************] 12.[**************]
6.[**************] 12.[**************]
7.[**************] 12.[**************]
8.[**************] 12.[**************]
9.[**************] 11.[**************]
10.[**************] 11.[**************]
11.[**************] 10.[**************]
12.[**************] 10.[**************]
13.[**************] 9.[**************]
14.[**************] 9.[**************]
15.[**************] 8.[**************]
16.[**************] 7.[**************]
17.[**************]
7.[**************]
18.[**************] 6.[**************]
19.[**************] 6.[**************]
20.[**************] 5.[**************]
21.[**************] 5.[**************]
22.[**************] 4.[**************]
23.[**************] 4.[**************]
24.[**************] 3.[**************]
25.[**************] 3.[**************]
26.[**************] 3.[**************]
27.[**************]
28.[**************]
29.[**************]
30.[**************]
31.[**************]
32.[**************]
33.[**************]
34.[**************]
35.[**************]
36.[**************]
37.[**************]
38.[**************]
39.[**************]
40.[**************]
41.[**************]
42.[**************]
43.[**************]
44.[**************]
45.[**************]
46.[**************]
47.[**************]
48.[**************]
49.[**************]
50.[**************]
51.[**************]
52.[**************]
53.[**************]
54.[**************]
55.[**************]
56.[**************]
57.[**************]
58.[**************]
○2输入命令 2.[**************] 2.[**************] 2.[**************] 2.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 1.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************] 0.[**************]
plot(t, vv,'g','LineWidth',2)
得燃料燃尽前火箭飞行速度变化曲线
运行
>> [t;vv]'
得火箭飞行中速度具体数据
ans =
0 0 1 12.[1**********]38 2 25.[1**********]97 3 37.[1**********]46 4 50.[1**********]84 5 63.[1**********]58 6 76.[1**********]68 7 88.[1**********]62 8 100.[1**********] 9 112.[1**********]8 10 123.[1**********]1 11 133.[1**********] 12
143.[1**********]3
13 153.[1**********]5 14 161.[1**********]2 15 170.[1**********]3 16 177.[1**********]6 17 184.[1**********]9 18 190.[1**********]2 19 196.[1**********]1 20 201.[1**********]5 21 206.[1**********]1 22 211.[1**********]1 23 214.[1**********]3 24 218.[1**********]7 25 221.[1**********]9 26 224.[1**********]6 27 227.[1**********]2 28 229.[1**********]4 29 231.[1**********]5 30 233.[1**********]9 31 235.[1**********]9 32 237.[1**********]1 33 239.[1**********]9 34 240.[1**********]7 35 242.[1**********]8 36 243.[1**********]4 37 244.[1**********]1 38 245.[1**********]5 39 247.[1**********]2 40 248.[1**********]4 41 249.[1**********]7 42 250.[1**********]8 43 251.[1**********]2 44 252.[1**********]3 45 253.[1**********]3 46 254.[1**********]2 47 255.[1**********]5 48 256.[1**********]4 49 257.[1**********]7 50 258.[1**********]6 51 259.[1**********]1 52 260.[1**********]9 53 261.[1**********] 54 262.[1**********]2 55 263.[1**********]7 56 264.[1**********]6
57 265.[1**********]2 58 266.[1**********]1 3输入命令 ○
>> plot(t, hh,'b','LineWidth',2)
得燃料燃尽前火箭飞行高度时间变化曲线
运行
>> [t;hh]'
得火箭飞行过程中高度具体数据
ans =
0 0 1 6.[1**********]771 2 25.[1**********]62 3 56.[1**********]49 4 100.[1**********]5 5 158.[1**********]4 6 227.[1**********]1 7 310.[1**********] 8
404.[1**********]7
10 628.[1**********]6 11 757.[1**********]4 12 895.[1**********]7 13 1044.[1**********] 14 1202.[1**********] 15 1368.[1**********] 16 1542.[1**********] 17 1723.1584828712 18 1910.[1**********] 19 2104.[1**********] 20 2304.[1**********] 21 2508.[1**********] 22 2717.[1**********] 23 2930.[1**********] 24 3147.[1**********] 25 3367.[1**********] 26 3590.[1**********] 27 3816.[1**********] 28 4045.[1**********] 29 4276.[1**********] 30 4509.[1**********] 31 4743.[1**********] 32 4980.[1**********] 33 5218.[1**********] 34 5458.[1**********] 35 5699.[1**********] 36 5942.[1**********] 37 6186.[1**********] 38 6431.[1**********] 39 6678.[1**********] 40 6926.[1**********] 41 7174.[1**********] 42 7424.[1**********] 43 7675.[1**********] 44 7928.[1**********] 45 8181.[1**********] 46 8435.508981705 47 8690.[1**********] 48 8947.[1**********] 49 9204.[1**********] 50 9462.[1**********] 51 9721.[1**********] 52 9982.2660536508
54 10505.8206883216 55 10769.0654797243 56 11033.2819028965 57 11298.4652041564 58 11564.6107884638