s=0; for i=1:n
if r(i)>180 s=s+1; end end pp=s/n
3、结果如下：pp =0.0486
实验六
6．1 用欧拉公式和四阶龙格-库塔法分别求解下列初值问题； 一、问题描述
6．1 用欧拉公式和四阶龙格-库塔法分别求解下列初值问题；
0.9y,y(0)?1;x?[0,1]1?2x
xy(2)y'??,y(0)?2;x?[0,1]21?x(1)y'?二、源程序及运行结果
(1)a.欧拉法
function euler(a,b,h,alpha) t=a:h:b; n=(b-a)/h; y(1)=alpha; for i=2:n+1
y(i)=y(i-1)+h*f(t(i-1),y(i-1));
sprintf('t=%3.1f,y=%9.7f',t(i),y(i)) end
①function r=f(t,y) r=(0.9*y)/(1+2*t);
在Mathlab Command Window里输入如下代码:euler(0,1,0.1,1) 得到：
ans =
t=0.1,y=1.0900000 ans =
t=0.2,y=1.1717500 ans =
t=0.3,y=1.2470768 ans =
t=0.4,y=1.3172249 ans =
t=0.5,y=1.3830861 ans =
t=0.6,y=1.4453250 ans =
t=0.7,y=1.5044519 ans =
t=0.8,y=1.5608688 ans =
t=0.9,y=1.6148989 ans =
t=1.0,y=1.6668064 ②
function r=f(x,y) r=-x*y/(1+x*x);
在Mathlab Command Window里输入如下代码 euler(0,1,0.1,2) 结果：
ans =
t=0.1,y=2.0000000 ans =
t=0.2,y=1.9801980 ans =
t=0.3,y=1.9421173 ans =
t=0.4,y=1.8886645 ans =
t=0.5,y=1.8235382 ans =
t=0.6,y=1.7505966 ans =
t=0.7,y=1.6733644 ans =
t=0.8,y=1.5947500 ans =
t=0.9,y=1.5169573 ans =
t=1.0,y=1.4415285 b. 四阶龙格-库塔法
fun=inline('0.9*y/(1+2*t)','t','y'); [t,y]=ode45(fun,[0,1],1); [t,y] ans =
0 1.0000
0.0250 1.0222 0.0500 1.0438
0.0750 0.1000 0.1250 1.0649 1.0855 1.1056
0.1500 1.1253 0.6000 1.4259 0.1750 1.1446 0.6250 1.4404 0.2000 1.1635 0.6500 1.4547 0.2250 1.1820 0.6750 1.4689 0.2500 1.2002 0.7000 1.4828 0.2750 1.2180 0.7250 1.4967 0.3000 1.2355 0.7500 1.5103 0.3250 1.2528 0.7750 1.5239 0.3500 1.2697 0.8000 1.5372 0.3750 1.2864 0.8250 1.5505 0.4000 1.3028 0.8500 1.5636 0.4250 1.3189 0.8750 1.5765 0.4500 1.3349 0.9000 1.5894 0.4750 1.3506 0.9250 1.6021 0.5000 1.3660 0.9500 1.6147 0.5250 1.3813 0.9750 1.6271 0.5500 1.3964 1.0000 1.6395 0.5750 1.4112
plot(t,y,'o-') %画图
（2）
odefun=inline('-x*y/(1+x*x)','x','y'); [x,y]=ode45(odefun,[0,1],1); [x,y]
ans =
0 1.0000
0.0250 0.9997 0.0500 0.9988 0.0750 0.9972 0.1000 0.9950 0.1250 0.9923 0.1500 0.9889 0.1750 0.9850 0.5250 0.8854 0.5500 0.8762 0.5750 0.8669 0.6000 0.8575 0.6250 0.8480 0.6500 0.8384 0.6750 0.8288 0.7000 0.8192 0.2000 0.9806 0.2250 0.9756 0.2500 0.9701 0.2750 0.9642 0.3000 0.9578 0.3250 0.9510 0.3500 0.9439 0.3750 0.9363 0.4000 0.9285 0.4250 0.9203 0.4500 0.9119 0.4750 0.9033 0.5000 0.8944
>> plot(x,y,'o-')
0.7250 0.7500 0.7750 0.8000 0.8250 0.8500 0.8750 0.9000 0.9250 0.9500 0.9750 1.0000 0.8096 0.8000 0.7904 0.7809 0.7714 0.7619 0.7526 0.7433 0.7341 0.7250 0.7160 0.7071
6．2 用求解常微分方程初值问题的方法计算积分上限函数 一、问题描述
6．2 用求解常微分方程初值问题的方法计算积分上限函数
y(x)??f(t)dt
ax的值，实际上是将上面表达式两端求导化为常微分方程形式，并用初值条件y(a)?0。
etdt的计算问题。 试用MATLAB中的指令ode23解决定积分?2?4t5二、源程序及运行结果
odefun=inline('2*pi*exp(x)/x','x','y');
ode23t(odefun,[4,5],0);