numerical analysis chapra 5ESM28.docx
- 文档编号:26382023
- 上传时间:2023-06-18
- 格式:DOCX
- 页数:56
- 大小:1.09MB
numerical analysis chapra 5ESM28.docx
《numerical analysis chapra 5ESM28.docx》由会员分享,可在线阅读,更多相关《numerical analysis chapra 5ESM28.docx(56页珍藏版)》请在冰豆网上搜索。
numericalanalysischapra5ESM28
CHAPTER28
28.1Thesolutionwiththe2nd-orderRK(Heunwithoutcorrector)canbelaidoutas
t
c
k1
c
k2
0
10
2
30
1
1.5
10
25
1.25
37.5
0.625
0.9375
20
34.375
0.78125
42.1875
0.390625
0.585938
30
40.23438
0.488281
45.11719
0.244141
0.366211
40
43.89648
0.305176
46.94824
0.152588
0.228882
50
46.1853
0.190735
48.09265
0.095367
0.143051
Forthe4th-orderRK,thesolutionis
t
c
k1
cmid
k2
cmid
k3
cend
k4
0
10
2
20
1.5
17.5
1.625
26.25
1.1875
1.572917
10
25.72917
1.213542
31.79688
0.910156
30.27995
0.986003
35.58919
0.72054
0.9544
20
35.27317
0.736342
38.95487
0.552256
38.03445
0.598278
41.25594
0.437203
0.579102
30
41.06419
0.446791
43.29814
0.335093
42.73965
0.363017
44.69436
0.265282
0.351382
40
44.57801
0.2711
45.93351
0.203325
45.59463
0.220268
46.78069
0.160965
0.213208
50
46.71009
0.164495
47.53257
0.123371
47.32695
0.133652
48.04662
0.097669
0.129369
Aplotofbothsolutionsalongwiththeanalyticalresultisdisplayedbelow:
28.2Themass-balanceequationscanbewrittenas
Selectedsolutionresults(Euler’smethod)aredisplayedbelow,alongwithaplotoftheresults.
t
c1
c2
c3
c4
c5
dc1/dt
dc2/dt
dc3/dt
dc4/dt
dc5/dt
0
0.0000
0.0000
0.0000
0.0000
0.0000
1.0000
0.0000
4.0000
0.0000
0.0000
1
1.0000
0.0000
4.0000
0.0000
0.0000
1.0800
0.2000
3.0000
0.3500
0.0400
2
2.0800
0.2000
7.0000
0.3500
0.0400
1.0872
0.3760
2.2600
0.5703
0.0848
3
3.1672
0.5760
9.2600
0.9203
0.1248
1.0488
0.5182
1.7138
0.6999
0.1307
4
4.2160
1.0942
10.9738
1.6201
0.2555
0.9839
0.6244
1.3113
0.7673
0.1752
5
5.1999
1.7186
12.2851
2.3874
0.4307
0.9051
0.6963
1.0147
0.7927
0.2165
76
13.2416
13.2395
18.6473
16.8515
12.9430
0.0002
0.0004
0.0001
0.0106
0.0179
77
13.2418
13.2399
18.6475
16.8621
12.9609
0.0002
0.0004
0.0001
0.0099
0.0168
78
13.2419
13.2403
18.6476
16.8720
12.9777
0.0001
0.0003
0.0001
0.0093
0.0158
79
13.2421
13.2406
18.6477
16.8813
12.9935
0.0001
0.0003
0.0001
0.0088
0.0149
80
13.2422
13.2409
18.6478
16.8901
13.0084
0.0001
0.0003
0.0001
0.0082
0.0140
Finally,MATLABcanbeusedtodeterminetheeigenvaluesandeigenvectors:
>>a=[.16-.06000;-.2.2000;0-.05.2500;00-.0875.125-.0375;-.04-.0200.06]
a=
0.1600-0.0600000
-0.20000.2000000
0-0.05000.250000
00-0.08750.1250-0.0375
-0.0400-0.0200000.0600
>>[v,d]=eig(a)
v=
000-0.10390.2604
000-0.1582-0.5701
00.81920-0.04360.6892
1.0000-0.57350.49970.4956-0.3635
000.86620.84660.0043
d=
0.12500000
00.2500000
000.060000
0000.06860
00000.2914
28.3Substitutingtheparametersintothedifferentialequationgives
Themid-pointmethodcanbeappliedwiththeresult:
Theresultsareapproachingasteady-statevalueof9.586267.
Challengequestion:
Thesteadystateform(i.e.,dc/dt=0)oftheequationis0=14.58333–0.08333c–0.15c2,whichcanbesolvedfor9.586267and–10.1418.Thus,thereisanegativeroot.
Ifweputintheinitialyvalueas10.1418(orhigherprecision),thesolutionwillstayatthenegativerootforawhileuntilroundofferrorsmayleadtothesolutiongoingunstable.Ifweuseaninitialvaluelessthanthenegativeroot,thesolutionwillalsogounstable.
However,ifwepickavaluethatisslightlyhigher(aspermachineprecision),itwillgravitatetowardsthepositiveroot.Forexampleifweuse10.14
Thiskindofbehaviorwillalsooccurforhigherinitialconditions.Forexample,usinganinitialconditionof16gives,
However,ifwestarttouseevenhigherinitialconditions,thesolutionwillagaingounstable.Forexample,ifweuseaninitialconditionof17,theresultis
Therefore,itlookslikeinitialconditionsroughlyintherangefrom-10.14uptoabout16.5willyieldstablesolutionsthatconvergeonthesteady-statesolutionof9.586.Notethatthisrangedependsonourchoiceofstepsize.
28.4Thefirststepsofthesolutionareshownbelowalongwithaplot.Noticethatavalueoftheinflowconcentrationattheendoftheinterval(cin-end)isrequiredtocalculatethek2’scorrectly.
t
cin
c
k1
cend
cin-end
k2
0
0
20
-1.2
17.6
8.534886
-0.54391
-0.87195
2
8.534886
18.25609
-0.58327
17.08955
15.24866
-0.11045
-0.34686
4
15.24866
17.56237
-0.13882
17.28472
20.52991
0.194711
0.027944
6
20.52991
17.61826
0.174699
17.96766
24.68428
0.402998
0.288848
8
24.68428
18.19595
0.3893
18.97455
27.95223
0.538661
0.46398
10
27.95223
19.12391
0.529699
20.18331
30.52289
0.620375
0.575037
Inthefollowingplot,theinflowconcentration(thinline)andtheoutflowconcentration(heavyline)arebothshown.
28.5Thesystemisasdepictedbelow:
(a)Themassofwaterinthetankcanbemodeledwithasimplemassbalance
WiththeinitialconditionthatV=1m3att=0,thisequationcanbeintegratedtoyield,
Thus,thetimetoemptythetank(M=0)canbecalculatedast=1/0.025=40hrs.
(b)Theconcentrationinthetankoverthisperiodcanbecomputedinseveralways.Thesimplestistocomputethemassofsaltinthetankovertimebysolvingthefollowingdifferentialequation:
(1)
wherem=themassofsaltinthetank.Thesaltconcentrationinthetank,s,istheratioofthemassofsalttothevolumeofwater
(2)
Thus,wecanintegrateEq.1todeterminethemassateachtimestepandthenuseEq.2todeterminethecorrespondingsaltconcentration.ThefirstfewstepsofthesolutionoftheODEusingEuler’smethodaretabulatedbelow.Inaddition,agraphoftheentiresolutionisalsodisplayed.Notethattheconcentrationasymptoticallyapproachesaconstantvalueof8,347.826g/m3asthesolutionevolves.
t
Vol
m
s
dm/dt
0
1
8000
8000
0
0.5
0.9875
8000
8101.266
-60.7595
1
0.975
7969.62
8173.969
-104.382
1.5
0.9625
7917.429
8225.901
-135.54
2
0.95
7849.659
8262.799
-157.679
2.5
0.9375
7770.819
8288.874
-173.324
Recognizethatasingularityoccursatt=40,becausethetankwouldbetotallyemptyatthispoint.
28.6Aheatbalanceforthespherecanbewrittenas
Theheatgaincanbetransformedintoavolumelossbyconsideringthelatentheatoffusion.
(1)
where=density0.917kg/m3andLf=latentheatoffusion333kJ/kg.Thevolumeandareaofaspherearecomputedby
(2)
Thesecanbecombinedwith
(1)toyield,
whereT=0oC.Thisequationcanbeintegratedalongwiththeinitialcondition,
toyieldtheresultingvolumeasafunctionoftime.Thisresultcanbeconvertedintodiameterusing
(2).Bothresultsareshownonthefollowingplotwhichindicatesthatittakesabout150sfortheicetomelt.
28.7Thesystemforthisproblemisstiff.Thus,theuseofasimpleexplicitRunge-Kuttaschemewouldinvolveusingaverysmalltimestepinordertomaintainastablesolution.Asolverdesignedforstiffsystemswasusedtogeneratethesolutionshownbelow.Twoviewsofthesolutionaregiven.Thefirstisfortheentiresolutiondomain.
Inaddition,wecanenlargetheinitialpartofthesolutiontoillustratethefasttransientsthatoccurasthesolutionmovesfromitsinitialconditionstoitsdominanttrajectories.
28.8Severalmethodscouldbeusedtoobtainasolutionforthisproblem(e.g.,finite-difference,shootingmethod,finite-element).Thefinite-differenceapproachisstraightforward:
Substitutingparametervaluesandcollectingtermsgives
Usingax=0.2cmthisequationcanbewrittenforalltheinteriornodes.TheresultinglinearsystemcanbesolvedwithanapproachliketheGauss-Seidelmethod.Thefollowingtableandgraphsummarizetheresults.
x
A
x
A
x
A
x
A
0
0.1
0.2
0.069544
1.2
0.011267
2.2
0.001779
3.2
0.000257
0.4
0.048359
1.4
0.007814
2.4
0.001224
3.4
0.000166
0.6
0.033621
1.6
0.005415
2.6
0.000840
3.6
9.93E-05
0.8
0.023368
1.8
0.003748
2.8
0.000574
3.8
4.65E-05
1
0.016235
2
0.002591
3
0.000389
4
0
28.9(a)ThetemperatureofthebodycanbedeterminedbyintegratingNewton’slawofcoolingtogive,
ThisequationcanbesolvedforK,
Substitutingthevaluesyields
Thetimeofdeathcanthenbecomputedas
Thus,thepersondied1.166hrspriortobeingdiscovered.Thenon-self-startingHeunyieldedthefollowingtimeseriesoftemperature.Forconvenience,wehaveredefinedthetimeofdeathast=0.:
28.10Theclassical4thorderRKmethodyields
28.11Theclassical4thorderRKmethodyields
t
C
Te
0
1
15
0.0625
0.941261
60.79675
0.125
0.885808
82.90298
0.1875
0.83356
92.37647
0.25
0.784373
95.1878
0.3125
0.738086
94.54859
0.375
0.694535
92.18087
0.4375
0.653562
89.00406
0.5
0.615016
85.50575
0.5625
0.578753
81.94143
0.625
0.544638
78.4421
0.6875
0.512542
75.07218
0.75
0.482346
71.86061
0.8125
0.453936
68.8176
0.875
0.427205
65.94362
0.9375
0.402055
63.23421
1
0.378391
60.68253
28.12InMATLAB,thefirststepistosetupafunctiontoholdthedifferentialequations:
functiondc=dcdtstiff(t,c)
dc=[-0.013*c
(1)-1000*c
(1)*c(3);-2500*c
(2)*c(3);-0.013*c
(1)-1000*c
(1)*c(3)-2500*c
(2)*c(3)];
Then,anODEsolverlikethefunctionode45canbeimplementedasin
>>tspan=[0,50];
>>y0=[1,1,0];
>>[t,y]=ode45(@dcdtstiff,tspan,y0);
Ifthisisdone,thesolutionwilltakearelativelylongtimetocomputetheresults.Incontrast,becauseitisexpresslydesignedtohandlestiffsystems,afunction
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- numerical analysis chapra 5ESM28 ESM28