dish（'Singularitylikely.'）
t
end
tout=tout（1:
k）;
yout=yout（1:
k,:
）;
流程图：
N
Y
N
Y
N
N
Y
Y
N
Y
N
Y
N
Y
N
Y
N
举例：
用欧拉法求y’=-y+x+1，y（0）=1。
解题思路：
首先建立例子所给函数的函数文件，调用函数myeuler，利用程序求解方程。
将欧拉法解和精确解比较，由方程f=-y+x+1可得到其精确解为y（x）=x+exp（-x）。
最后在同一图窗中分别画出两图。
程序代码：
f．m
functionf=f（x,y）
f=-y+x+1;
>>[x,y]=myeuler（'f',0,1,1）;%利用程序求解方程
>>y1=x+exp（-x）;%方程f=-y+x+1的精确解
>>plot（x,y,'-b',x,y1,'-r'）%在同一图窗将欧拉法解和精确解的图画出
>>legend（'欧拉法','精确解'）
例题流程图：
2、编程解决以下科学计算问题
（1）
解题思路：
◆建模：
材料力学中从弯矩求转角要经过一次不定积分,而从转角求挠度又要经过一次不定积分,在MATLAB中这却是非常简单的问题.可用cumsum函数作近似的不定积分,还可用更精确的函数cumtrapz来做不定积分。
本题用cumsum函数来做.解题的关键还是在于正确地列写弯矩方程。
本例中弯矩为
程序代码：
>>L=1;P=1000;L1=1;%给出常数
E=200*10^9;I=2*10^-5;
x=linspace（0,L,101）;dx=L/100;%将x分100段
n1=L1/dx+1;%确定x=L1处对应的下标
M1=-P*（L1-x（1:
n1））;%第一段弯矩赋值
M2=zeros（1,101-n1）;%第二段弯矩赋值（全为零）
M=[M1,M2];%全梁的弯矩
A=cumsum（M）*dx/（E*I）%对弯矩积分求转角
Y=cumsum（A）*dx%对转角积分求挠度
A=
1.0e-003*
Columns1through9
-0.0025-0.0050-0.0074-0.0098-0.0122-0.0146-0.0170-0.0193-0.0216
Columns10through18
-0.0239-0.0261-0.0283-0.0305-0.0327-0.0349-0.0370-0.0391-0.0412
Columns19through27
-0.0432-0.0452-0.0472-0.0492-0.0512-0.0531-0.0550-0.0569-0.0587
Columns28through36
-0.0605-0.0623-0.0641-0.0659-0.0676-0.0693-0.0710-0.0726-0.0742
Columns37through45
-0.0759-0.0774-0.0790-0.0805-0.0820-0.0835-0.0849-0.0863-0.0877
Columns46through54
-0.0891-0.0905-0.0918-0.0931-0.0944-0.0956-0.0968-0.0980-0.0992
Columns55through63
-0.1004-0.1015-0.1026-0.1037-0.1047-0.1057-0.1067-0.1077-0.1087
Columns64through72
-0.1096-0.1105-0.1114-0.1122-0.1130-0.1138-0.1146-0.1154-0.1161
Columns73through81
-0.1168-0.1175-0.1181-0.1187-0.1194-0.1199-0.1205-0.1210-0.1215
Columns82through90
-0.1220-0.1224-0.1229-0.1232-0.1236-0.1240-0.1243-0.1246-0.1249
Columns91through99
-0.1251-0.1253-0.1255-0.1257-0.1259-0.1260-0.1261-0.1262-0.1262
Columns100through101
-0.1262-0.1262
Y=
1.0e-004*
Columns1through9
-0.0002-0.0007-0.0015-0.0025-0.0037-0.0052-0.0069-0.0088-0.0109
Columns10through18
-0.0133-0.0159-0.0188-0.0218-0.0251-0.0286-0.0323-0.0362-0.0403
Columns19through27
-0.0446-0.0492-0.0539-0.0588-0.0639-0.0692-0.0747-0.0804-0.0863
Columns28through36
-0.0924-0.0986-0.1050-0.1116-0.1184-0.1253-0.1324-0.1397-0.1471
Columns37through45
-0.1547-0.1624-0.1703-0.1783-0.1865-0.1949-0.2034-0.2120-0.2208
Columns46through54
-0.2297-0.2388-0.2479-0.2572-0.2667-0.2762-0.2859-0.2957-0.3057
Columns55through63
-0.3157-0.3258-0.3361-0.3465-0.3569-0.3675-0.3782-0.3890-0.3998
Columns64through72
-0.4108-0.4218-0.4330-0.4442-0.4555-0.4669-0.4784-0.4899-0.5015
Columns73through81
-0.5132-0.5249-0.5367-0.5486-0.5606-0.5726-0.5846-0.5967-0.6088
Columns82through90
-0.6210-0.6333-0.6456-0.6579-0.6703-0.6827-0.6951-0.7075-0.7200
Columns91through99
-0.7325-0.7451-0.7576-0.7702-0.7828-0.7954-0.8080-0.8206-0.8332
Columns100through101
-0.8459-0.8585
>>subplot（3,1,1）,plot（x,M）,grid%绘弯矩图
subplot（3,1,2）,plot（x,A）,grid%绘弯矩图
subplot（3,1,3）,plot（x,Y）,grid%绘弯矩图
流程图
（2）计算积分
①
解题思路：
exp（-x^2）是不可积函数，matlab中int积分显示不出积分结果，用到vpa函数控制其结果精度，表示出来。
程序：
>>symsx
>>t=vpa（int（exp（-x.^2）./（1+x.^2）,-inf,+inf）,5）
结果：
t=1.3433
②
解题思路：
先建立内置函数，然后调用quad函数求积分。
程序：
>>fun=@（x）tan（x）./x.^（0.7）;
>>quad（'fun',0,1）
结果：
ans=0.9063
③
解题思路：
先建立内置函数，然后调用quad函数求积分。
程序：
>>fun=inline（'exp（x）./（1-x.^2）.^0.5'）;
>>quad（'fun',0,1）
结果：
ans=3.1044
④
解题思路：
这是一个二重积分，一般的方法是把二重积分化为二次积分，但由于二次积分命令
int（int（1+x+y.^2,y,-sqrt（2*x-x.^2）,sqrt（2*x-x.^2））,x,0,1）
无法求出结果，故将其转化成极坐标。
积分函数为
r.*（1+r.*cos（a）+r.^2.*（sin（a））.^2）
其中a的范围：
-pi/2到pi/2；r的范围：
0到2.*cos（a）。
程序：
>>symsar
>>int（int（r+r.^2.*cos（a）+r.^3.*（sin（a））.^2,r,0,2.*cos（a））,a,-pi/2,pi/2）
结果：
ans=（9*pi）/4
友情提示：
方案范本是经验性极强的领域，本范文无法思考和涵盖全面，供参考!
最好找专业人士起草或审核后使用。
升级会员 