数学建模实验.docx
《数学建模实验.docx》由会员分享，可在线阅读，更多相关《数学建模实验.docx（10页珍藏版）》请在冰点文库上搜索。
数学建模实验
数学建模课程实验报告
专题实验7
班级
数财系1班
学号
2011040123
丛文
实验题目
常微分方程数值解
实验目的
1．掌握用MATLAB求微分方程初值问题数值解的方法；
2．通过实例学习微分方程模型解决简化的实际问题；
3．了解欧拉方法和龙格库塔方法的基本思想。
实验容
（包括分析过程、方法、和代码，结果）
1.用欧拉方法和龙格库塔方法求下列微分方程初值问题的数值解，画出解的图形，对结果进行分析比较
解；M文件
functionf=f（x,y）
f=y+2*x;
程序；
clc;clear;
a=0;b=1;%求解区间
[x1,y_r]=ode45（'f',[ab],1）;%调用龙格库塔求解函数求解数值解；
%%以下利用Euler方法求解
y
（1）=1;N=100;h=（b-a）/N;
x=a:
h:
b;
fori=1:
N
y（i+1）=y（i）+h*f（x（i）,y（i））;
end
figure
（1）
plot（x1,y_r,'r*',x,y,'b+',x,3*exp（x）-2*x-2,'k-'）;%数值解与真解图
title（'数值解与真解图'）;
legend（'RK4','Euler','真解'）;
xlabel（'x'）;ylabel（'y'）;
figure
（2）
plot（x1,abs（y_r-（3*exp（x1）-2*x1-2））,'k-'）;%龙格库塔方法的误差
title（'龙格库塔方法的误差'）
xlabel（'x'）;ylabel（'Error'）;
figure（3）
plot（x,abs（y-（3*exp（x）-2*x-2））,'r-'）%Euler方法的误差
title（'Euler方法的误差'）
xlabel（'x'）;ylabel（'Error'）;
4.单摆运动是一个我们熟悉的物理模型，可以看作工程技术中一些振动问题的简化，图8中一根长l的（无弹性的）细线，一端固定，另一端悬挂一质量为m的小球，在重力作用下小球处于竖直的平衡位置，使小球偏离平衡位置一根小的角度
，然后让它无初速度的放开，小球就会沿圆弧摆动，在不考虑空气阻力的情况下建立
关于时间t的微分方程，设l=25cm，在
等于
和
两种情况下求方程的数值解，并与近似解
比较。
解；实验原理与数学模型：
在小球摆动过程中的任一位置θ，小球所受重力沿运动轨迹方向的分力为-mgsinθ（负号表示力的方向与θ的正方向相反），利用牛顿第二定律即得微分方程
描述单摆运动规律的微分方程
（1）是2阶微分方程，无解析解，但可用Matlab或其它软件编程求其数值解，但都需要先将它化成方程组的形式。
令
则微分方程
（1）化为
初始条件转化为
在前面的两式中，g=9.8，l=0.25，x10为10o=0.1745（弧度）及30o=0.5236（弧度）两种情况.
周期
根据上原理，可以建立模型，用matlab编程求解，过程如下，程序以单摆的两个周期来计算作图。
M文件；
functiondx=danbai（t,x）
g=9.8;l=0.25;
dx=[x
（2）;-g/l*sin（x
（1））];
程序；
一.当
时
ts=0:
0.05:
2;
a0=0.1745;
x0=[a0,0];
[t,x]=ode23（danbai,ts,x0）;
y=a0*cos（sqrt（40）.*t）;
[t,x（:
1）,y]
subplot（1,2,2）,plot（t,x（:
1）,'-k*'）,title（'摆角10度数值解'）
subplot（1,2,1）,plot（t,y,'b*'）,plot（t,y,'-r*'）,title（'摆角10度近似解'）
ans=
00.17450.1745
0.05000.16610.1658
0.10000.14150.1407
0.15000.10330.1017
0.20000.05500.0525
0.25000.0015-0.0018
0.3000-0.0522-0.0560
0.3500-0.1009-0.1046
0.4000-0.1397-0.1429
0.4500-0.1649-0.1669
0.5000-0.1742-0.1745
0.5500-0.1667-0.1647
0.6000-0.1430-0.1386
0.6500-0.1054-0.0987
0.7000-0.0577-0.0491
0.7500-0.00440.0054
0.80000.04940.0594
0.85000.09830.1075
0.90000.13770.1449
0.95000.16370.1679
1.00000.17380.1744
二.当
时，
ts=0:
0.05:
2;
a0=0.5236;
x0=[a0,0];
[t,x]=ode23（danbai,ts,x0）;
y=a0*cos（sqrt（40）.*t）;
[t,x（:
1）,y]
subplot（1,2,2）,plot（t,x（:
1）,'-k*'）,title（'摆角30度数值解'）
subplot（1,2,1）,plot（t,y,'b*'）,plot（t,y,'-r*'）,title（'摆角30度近似解'）
ans=
00.52360.5236
0.05000.49930.4976
0.10000.42830.4223
0.15000.31700.3051
0.20000.17530.1577
0.25000.0169-0.0054
0.3000-0.1431-0.1680
0.3500-0.2893-0.3139
0.4000-0.4076-0.4286
0.4500-0.4874-0.5009
0.5000-0.5217-0.5235
0.5500-0.5075-0.4942
0.6000-0.4460-0.4158
0.6500-0.3425-0.2963
0.7000-0.2064-0.1473
0.7500-0.05050.0162
0.80000.11020.1782
0.85000.26020.3225
0.90000.38520.4348
0.95000.47360.5039
1.00000.51770.5232
1.05000.51360.4905
1.10000.46170.4092
1.15000.36650.2873
1.20000.23650.1369
1.25000.0838-0.0271
1.3000-0.0770-0.1883
1.3500-0.2303-0.3309
1.4000-0.3613-0.4407
1.4500-0.4580-0.5068
1.5000-0.5116-0.5226
1.5500-0.5177-0.4866
1.6000-0.4755-0.4023
1.6500-0.3889-0.2782
1.7000-0.2654-0.1264
1.7500-0.11640.0379
1.80000.04370.1984
1.85000.19950.3393
1.90000.33610.4465
1.95000.44050.5094
2.00000.50360.5218
从数据可以看出，角度为10o时精确（数值）解与近似解相差不大，而初始角度为30o时，随着时间的增加差别很大
实验结果分析或者实验总结和体会
龙格库塔方法和Euler方法求解常微分方程都能获得比较好的数值解，相比较而言龙格库塔方法的数值解的精度远远要比Euler方法的数值解的精度高。
升级会员 