实验四 各类方程的求解matlab实验报告.docx
- 文档编号:12717820
- 上传时间:2023-04-21
- 格式:DOCX
- 页数:53
- 大小:280.40KB
实验四 各类方程的求解matlab实验报告.docx
《实验四 各类方程的求解matlab实验报告.docx》由会员分享,可在线阅读,更多相关《实验四 各类方程的求解matlab实验报告.docx(53页珍藏版)》请在冰豆网上搜索。
实验四各类方程的求解matlab实验报告
实验四各类方程的求解
姓名:
学号:
10000实验日期:
2013.4.17
实验目的:
掌握求解各类方程的MATLAB命令
实验项目:
(1)非线性方程求根;非线性方程求极值;
(2)数值积分;
(3)常微分方程求解。
实验背景:
客观世界中众多现象可以划归为各类方程,特别是微分方程的求解问题;因而如何利用MATLAB求解各类方程的根、极值以及解是一类重要的课题。
实验具体过程:
1、题目:
(方程求根)
(i)求方程
的正根;
(ii)Newton迭代法是一种速度很快的迭代方法,但是它需要预先求得导函数,若用差商代替导数,可得下列弦截法
。
这一迭代法需要两个初值
,编写一个通用的弦截法计算机程序并用以解(i)。
解题思路:
(i)由于该方程是非线性、非多项式的方程,因而求解只能利用fzero或者fsolve命令;
(ii)和Newton迭代法不同,弦截法涉及两次递归,求解过程应避免出现
很小的情形。
对实验题目的解答:
(i)>>fun=inline('x.*log(sqrt(x.^2-1)+x)-sqrt(x.^2-1)-0.5*x');
>>fplot(fun,[1,10])
由此,可以看到零点大约在x=1附近。
下面使用fsolve命令进行计算精确的零点位置。
>>[x,f,h]=fsolve(fun,1)
Optimizerappearstobeconvergingtoapointwhichisnotaroot.
NormofrelativechangeinXislessthanmax(options.TolX^2,eps)but
sum-of-squaresoffunctionvaluesisgreaterthanorequaltosqrt(options.TolFun)
Tryagainwithanewstartingguess.
x=
0.7500-0.0000i
f=
-0.3750-0.1194i
h=
-2
>>[x,f,h]=fsolve(fun,1.2)
Optimizationterminated:
first-orderoptimalityislessthanoptions.TolFun.
x=
2.1155
f=
1.0709e-006
h=
1
(ii)%弦截法tps.m
functiony=tps(fun,x0,x1)
e=1e-6;
ifabs(x1-x0)<2*e,x1=x0+2*e;end
y=x1;
whileabs(x1-x0)>e
y=x1-(x1-x0)/(feval(fun,x1)-feval(fun,x0))*feval(fun,x1);
x0=x1;
x1=y;
end
>>tps(fun,1.1,1.2)
ans=
2.1155
改进或思考:
由(i)的运算结果可以看到,对于非线性问题,初始迭代点的选取非常重要。
一旦初始迭代数据过于远离精确值,将导致错误结果。
2、题目:
(函数的的极值问题)
考虑函数
(i)作出
在
的图,观察极值点的位置;
(ii)用MATLAB函数fminsearch求极值点和极值。
对实验题目的解答:
(i)>>clear;
>>xa=-2:
0.1:
1;ya=-7:
0.1:
1;
>>[x,y]=meshgrid(xa,ya);
>>z=y.^3/9+3*x.^2.*y+9*x.^2+y.^2+x.*y+9;
>>mesh(x,y,z)
从图形,我们可以大致估计函数有三个极大值,分别对应(-2,15),(0,22),(5,40),三个极小值,分别对应(-6,-20),(2,5),(-7,0)。
(ii)
下面利用fminsearch确定函数极值得精确位置。
极大值:
>>fun='x
(2)^3/9+3*x
(1)^2*x
(2)+9*x
(1)^2+x
(2)^2+x
(1)*x
(2)+9';
>>[X1,Y1]=fminsearch(['-',fun],[-215]);
[X2,Y2]=fminsearch(['-',fun],[022]);
[X3,Y3]=fminsearch(['-',fun],[540]);
Exiting:
Maximumnumberoffunctionevaluationshasbeenexceeded
-increaseMaxFunEvalsoption.
Currentfunctionvalue:
-130445905685433800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
Exiting:
Maximumnumberoffunctionevaluationshasbeenexceeded
-increaseMaxFunEvalsoption.
Currentfunctionvalue:
-1267583789692283300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
Exiting:
Maximumnumberoffunctionevaluationshasbeenexceeded
-increaseMaxFunEvalsoption.
Currentfunctionvalue:
-3096693623978009000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
>>[X1,-Y1],[X2,-Y2],[X3,-Y3]
ans=
1.0e+134*
0.00000.00001.3045
ans=
1.0e+135*
-0.00000.00001.2676
ans=
1.0e+135*
-0.00000.00003.0967
极小值:
>>[X1,Z1]=fminsearch(fun,[-6-20]);
[X2,Z2]=fminsearch(fun,[25]);
[X3,Z3]=fminsearch(fun,[-70]);
Exiting:
Maximumnumberoffunctionevaluationshasbeenexceeded
-increaseMaxFunEvalsoption.
Currentfunctionvalue:
-1032540430221527300000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000.000000
>>[X1,Z1],[X2,Z2],[X3,Z3]
ans=
1.0e+135*
-0.0000-0.0000-1.0325
ans=
-0.00000.00009.0000
ans=
-0.00010.00109.0000
改进或思考:
本题之所以运行结果与图形差异较大,是因为没有限制自变量的范围。
想得到精确地解,须对原函数进行补充定义,如考察
的极小(大)值。
3、题目:
(Simpson积分法)编制一个定步长Simpson法数值积分程序。
计算公式为
,
其中
为偶数,
,
,并取
,求解积分
的数值解。
解题思路:
定义一M函数文件来描述Simpson积分法,这里可选择步长h,步数n,精度e作为函数的参数。
对实验题目的解答:
functionsimpson(fun,a,b,m,e)
disp('请确认您依次输入的是函数,左端点,右端点,步数,精度。
')
n=2*m;h=(b-a)/n;S=h*(feval(fun,a)+feval(fun,b))/3;
fori=2:
2:
n
S=S+4*feval(fun,a+(i-1)*h);
end
fori=3:
2:
n-1
S=S+2*feval(fun,a+(i-1)*h);
end
S
>>fun=inline('exp(-x.^2/2)/sqrt(2*pi)');
>>simpson(fun,0,1,5,1e-6)
请确认您依次输入的是函数,左端点,右端点,步数,精度。
S=
9.6208
3、题目:
(电视机价格)由于市场竞争的影响,电视机售价p越高,销售量x就会越低,
,其中M为最大需求量,a为价格系数。
另一方面销售量越大,每台电视机成本c就会越低,
,其中
是只生产一台电视机时的成本,k为规模系数。
应如何确定电视机售价才能获得最大利润?
解题思路:
建立电视机总理论的函数表达式y=f(p;M,a,c0,k);通过求f’(p;M,a,c0,k)=0得到最大利润对应的售价。
本题的解题思路有两种方案:
1、利用定义M函数文件确认不同参数情况下售价的取法;2、利用符号运算直接求得售价的取法。
对实验题目的解答:
设电视机的总理论为y=f(p;M,a,c0,k);根据题意,
。
为此,我们需求出
的实根即可。
利用第七章所学的符号运算,我们可以求出
。
>>symspMac0kP;
>>P=M*exp(-a*p)*(p-c0+k*log(M*exp(-a*p)));
>>Q=diff(P,p)
Q=
-M*a*exp(-a*p)*(p-c0+k*log(M*exp(-a*p)))+M*exp(-a*p)*(1-k*a)
下面我们利用两种方式求解该问题:
(1)编写一M函数文件
%电视机价格问题price.m
functionq=price(M,a,c0,k,Tp)
P=@(p)M*exp(-a*p)*(p-c0+k*log(M*exp(-a*p)));
Q=@(p)-M*a*exp(-a*p).*(p-c0+k*log(M*exp(-a*p)))+M*exp(-a*p).*(1-k*a);
fplot(Q,[c0,10*c0]);
holdon;
fplot(P,[c0,10*c0])
[p,Q0,h]=fsolve(Q,Tp)
feval(P,p)
(2)符号运算
>>q=solve(Q,p)
q=
-(a*c0+k*a*log(1/M)+1-k*a)/a/(-1+k*a)
>>subs(P,p,q)%计算最大利润
ans=
M*exp((a*c0+k*a*log(1/M)+1-k*a)/(-1+k*a))*(-(a*c0+k*a*log(1/M)+1-k*a)/a/(-1+k*a)-c0+k*log(M*exp((a*c0+k*a*log(1/M)+1-k*a)/(-1+k*a))))
现对[M,a,c0,k]=[1,0.5,1,1],分别利用上述结果进行计算。
(1)>>price(1,0.5,1,1,1)
Optimizationterminated:
first-orderoptimalityislessthanoptions.TolFun.
p=
4.0000
Q0=
2.0734e-008
h=
1
ans=
0.1353
>>price(1,0.5,1,1,4)
Optimizationterminated:
first-orderoptimalityislessthanoptions.TolFun.
p=
4
Q0=
0
h=
1
ans=
0.1353
(2)>>M=1;a=0.5;c0=1;k=1;
>>-(a*c0+k*a*log(1/M)+1-k*a)/a/(-1+k*a)
ans=
4
>>M*exp((a*c0+k*a*log(1/M)+1-k*a)/(-1+k*a))*(-(a*c0+k*a*log(1/M)+1-k*a)/a/(-1+k*a)-c0+k*log(M*exp((a*c0+k*a*log(1/M)+1-k*a)/(-1+k*a))))
ans=
0.1353
由此可知,两者算出的结果是一致的。
改进或思考:
对于编写者而言,显然利用符号运算更为简便;然而对于过于复杂的表达式,符号计算可能失效。
因而两种方式各有利弊。
4、题目:
(刚性方程组求解)求解刚性方程组。
。
解题思路:
由于该方程组中一次项系数远大于零此项系数,该方程组为一刚性方程组。
可开用ode15s命令求解。
对实验题目的解答:
%ex4.m
functionf=ex4(t,y)
f
(1)=-1000.25*y
(1)+999.75*y
(2)+0.5;
f
(2)=999.75*y
(1)-1000.25*y
(2)+0.5;
f=f(:
);
>>[t,y]=ode15s(@ex4,[0,50],[1,-1]);plot(t,y)
5、题目:
(肿瘤生长)肿瘤大小V生长的速率与V的a次方成正比,其中a为形状参数,
;而其比例系数K随时间减小,减小速率又与当时的K值成正比,比例系数为环境参数b。
设某肿瘤参数a=1,b=0.1,K的初始值为2,V的初始值为1。
问
(i)此肿瘤生长不会超过多大?
(ii)过多长时间肿瘤大小翻一倍?
(iii)何时肿瘤生长速率由递增转为递减?
(iv)若参数a=2/3呢?
解题思路:
建立关于肿瘤大小V和比例系数K的常微分方程组,这里可以利用常微分方程组的数值解法(第六章)求出相应问题的近似结果;或者直接利用第七章符号运算的方法求出相应问题的精确解。
对实验题目的解答:
首先建立相应的常微分方程组
(1)符号运算
>>s=dsolve('DV=K*V','DK=-0.1*K','V(0)=1','K(0)=2')
s=
K:
[1x1sym]
V:
[1x1sym]
>>s.K,s.V
ans=
2*exp(-1/10*t)
ans=
exp
(2)^10*exp(-20*exp(-1/10*t))
由此,知道肿瘤大小不会超过exp
(2)^10。
>>symst;
>>T0=solve(exp
(2)^10*exp(-20*exp(-1/10*t))-2,t)
T0=
-10*log(-1/20*log(16777216/4069860639536131))
当t=T0时,肿瘤大小增大一倍。
>>u=diff(exp
(2)^10*exp(-20*exp(-1/10*t)),t,2)
u=
-4069860639536131/41943040*exp(-1/10*t)*exp(-20*exp(-1/10*t))+4069860639536131/2097152*exp(-1/10*t)^2*exp(-20*exp(-1/10*t))
>>T1=solve(u,t)
T1=
10*log(20)
当T=T1时,肿瘤生长速率由递增变为递减
注:
T0,T1用十进制数表示为
>>t0=vpa(T0),t1=vpa(T1)
t0=
.35272172300235666784997847812874
t1=
29.957322735539909934352235761425
将上述过程用a=2/3再次进行计算发现,符号运算进入长时间计算状态,失败!
(2)利用常微分方程组数值解的方法求解。
%肿瘤问题cancer.m
functionf=cancer(t,x)
f
(1)=x
(1)*x
(2);
f
(2)=-0.1*x
(2);
f=f(:
);
我们发现V和K的数量级悬殊太大,导致无法区分K的变化,为此我们重新单独输出V和K。
>>[t,y]=ode45(@cancer,[050],[1;2]);
>>t
t=
0
0.0251
0.0502
0.0754
0.1005
0.2121
0.3237
0.4354
0.5470
0.6857
0.8244
0.9631
1.1018
1.2544
1.4071
1.5597
1.7124
1.8761
2.0398
2.2035
2.3672
2.5420
2.7168
2.8915
3.0663
3.2533
3.4402
3.6272
3.8141
4.0149
4.2156
4.4163
4.6170
4.8335
5.0500
5.2666
5.4831
5.7180
5.9529
6.1877
6.4226
6.6791
6.9355
7.1920
7.4485
7.7306
8.0127
8.2948
8.5770
8.8901
9.2033
9.5164
9.8295
10.1808
10.5320
10.8833
11.2345
11.6335
12.0325
12.4316
12.8306
13.2910
13.7514
14.2118
14.6723
15.2139
15.7554
16.2970
16.8386
17.4912
18.1437
18.7963
19.4488
20.2586
21.0684
21.8782
22.6879
23.7258
24.7636
25.8014
26.8392
28.0892
29.3392
30.5892
31.8392
33.0892
34.3392
35.5892
36.8392
38.0892
39.3392
40.5892
41.8392
43.0892
44.3392
45.5892
46.8392
47.6294
48.4196
49.2098
50.0000
>>[t,y]=ode45(@cancer,[01000],[1;2])
t=
1.0e+003*
0
0.0000
0.0001
0.0001
0.0001
0.0002
0.0003
0.0004
0.0005
0.0007
0.0008
0.0010
0.0011
0.0013
0.0014
0.0016
0.0017
0.0019
0.0020
0.0022
0.0024
0.0025
0.0027
0.0029
0.0031
0.0033
0.0034
0.0036
0.0038
0.0040
0.0042
0.0044
0.0046
0.0048
0.0051
0.0053
0.0055
0.0057
0.0060
0.0062
0.0064
0.0067
0.0069
0.0072
0.0074
0.0077
0.0080
0.0083
0.0086
0.0089
0.0092
0.0095
0.0098
0.0102
0.0105
0.0109
0.0112
0.0116
0.0120
0.0124
0.0128
0.0133
0.0138
0.0142
0.0147
0.0152
0.0158
0.0163
0.0168
0.0175
0.0181
0.0188
0.0194
0.0203
0.0211
0.0219
0.0227
0.0237
0.0248
0.0258
0.0268
0.0282
0.0295
0.0309
0.0323
0.0340
0.0357
0.0374
0.0390
0.0410
0.0430
0.0450
0.0470
0.0489
0.0509
0.0529
0.0548
0.0568
0.0588
0.0607
0.0627
0.0647
0.0666
0.0686
0.0706
0.0725
0.0745
0.0765
0.0784
0.0804
0.0824
0.0843
0.0863
0.0884
0.0904
0.0925
0.0945
0.0969
0.0994
0.1018
0.1042
0.1070
0.1098
0.1126
0.1154
0.1188
0.1222
0.1255
0.1289
0.1331
0.1372
0.1414
0.1455
0.1508
0.1562
0.1615
0.1668
0.1739
0.1810
0.1882
0.1953
0.2050
0.2146
0.2243
0.2340
0.2449
0.2558
0.2667
0.2776
0.2840
0.2904
0.2969
0.3033
0.3097
0.3162
0.3226
0.3291
0.3378
0.3466
0.3554
0.3642
0.3751
0.3860
0.3968
0.4077
0.4147
0.4217
0.4288
0.4358
0.4428
0.4498
0.4568
0.4638
0.4727
0.4817
0.4906
0.4995
0.5098
0.5201
0.5305
0.5408
0.5479
0.5550
0.5620
0.5691
0.5762
0.5833
0.5903
0.5974
0.6062
0.6150
0.6239
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 实验四 各类方程的求解matlab实验报告 实验 各类 方程 求解 matlab 报告