Two-body decays of Lambda_b

a r X i v :h e p -e x /0207007v 1 1 J u l 2002BELLEBelle Prerpint 2002-18KEK Preprint 2002-59Study of B →ρπdecays at BelleBelle Collaboration A.Gordon u ,Y.Chao z ,K.Abe h ,K.Abe aq ,N.Abe at ,R.Abe ac ,T.Abe ar ,Byoung Sup Ahn o ,H.Aihara as ,M.Akatsu v ,Y.Asano ay ,T.Aso aw ,V.Aulchenko b ,T.Aushev ℓ,A.M.Bakich an ,Y.Ban ag ,A.Bay r ,I.Bedny b ,P.K.Behera az ,jak m ,A.Bondar b ,A.Bozek aa ,M.Braˇc ko t ,m ,T.E.Browder g ,B.C.K.Casey g ,M.-C.Chang z ,P.Chang z ,B.G.Cheon am ,R.Chistov ℓ,Y.Choi am ,Y.K.Choi am ,M.Danilov ℓ,L.Y.Dong j ,J.Dragic u ,A.Drutskoy ℓ,S.Eidelman b ,V.Eiges ℓ,Y.Enari v ,C.W.Everton u ,F.Fang g ,H.Fujii h ,C.Fukunaga au ,N.Gabyshev h ,A.Garmash b ,h ,T.Gershon h ,B.Golob s ,m ,R.Guo x ,J.Haba h ,T.Hara ae ,Y.Harada ac ,N.C.Hastings u ,H.Hayashii w ,M.Hazumi h ,E.M.Heenan u ,I.Higuchi ar ,T.Higuchi as ,L.Hinz r ,T.Hokuue v ,Y.Hoshi aq ,S.R.Hou z ,W.-S.Hou z ,S.-C.Hsu z ,H.-C.Huang z ,T.Igaki v ,Y.Igarashi h ,T.Iijima v ,K.Inami v ,A.Ishikawa v ,H.Ishino at ,R.Itoh h ,H.Iwasaki h ,Y.Iwasaki h ,H.K.Jang a ℓ,J.H.Kang bc ,J.S.Kang o ,N.Katayama h ,Y.Kawakami v ,N.Kawamura a ,T.Kawasaki ac ,H.Kichimi h ,D.W.Kim am ,Heejong Kim bc ,H.J.Kim bc ,H.O.Kim am ,Hyunwoo Kim o ,S.K.Kim a ℓ,T.H.Kim bc ,K.Kinoshita e ,S.Korpar t ,m ,P.Krokovny b ,R.Kulasiri e ,S.Kumar af ,A.Kuzmin b ,Y.-J.Kwon bc ,nge f ,ai ,G.Leder k ,S.H.Lee a ℓ,J.Li ak ,A.Limosani u ,D.Liventsevℓ,R.-S.Lu z,J.MacNaughton k,G.Majumder ao, F.Mandl k,D.Marlow ah,S.Matsumoto d,T.Matsumoto au,W.Mitaroffk,K.Miyabayashi w,Y.Miyabayashi v,H.Miyake ae,H.Miyata ac,G.R.Moloney u,T.Mori d,T.Nagamine ar,Y.Nagasaka i,T.Nakadaira as,E.Nakano ad, M.Nakao h,J.W.Nam am,Z.Natkaniec aa,K.Neichi aq, S.Nishida p,O.Nitoh av,S.Noguchi w,T.Nozaki h,S.Ogawa ap, T.Ohshima v,T.Okabe v,S.Okuno n,S.L.Olsen g,Y.Onuki ac, W.Ostrowicz aa,H.Ozaki h,P.Pakhlovℓ,H.Palka aa,C.W.Park o,H.Park q,L.S.Peak an,J.-P.Perroud r, M.Peters g,L.E.Piilonen ba,J.L.Rodriguez g,F.J.Ronga r, N.Root b,M.Rozanska aa,K.Rybicki aa,H.Sagawa h,S.Saitoh h,Y.Sakai h,M.Satapathy az,A.Satpathy h,e,O.Schneider r,S.Schrenk e,C.Schwanda h,k,S.Semenovℓ,K.Senyo v,R.Seuster g,M.E.Sevior u,H.Shibuya ap,V.Sidorov b,J.B.Singh af,S.Staniˇc ay,1,M.Stariˇc m,A.Sugi v, A.Sugiyama v,K.Sumisawa h,T.Sumiyoshi au,K.Suzuki h,S.Suzuki bb,S.Y.Suzuki h,T.Takahashi ad,F.Takasaki h, K.Tamai h,N.Tamura ac,J.Tanaka as,M.Tanaka h,G.N.Taylor u,Y.Teramoto ad,S.Tokuda v,S.N.Tovey u,T.Tsuboyama h,T.Tsukamoto h,S.Uehara h,K.Ueno z, Y.Unno c,S.Uno h,hiroda h,G.Varner g,K.E.Varvell an,C.C.Wang z,C.H.Wang y,J.G.Wang ba,M.-Z.Wang z,Y.Watanabe at,E.Won o,B.D.Yabsley ba,Y.Yamada h, A.Yamaguchi ar,Y.Yamashita ab,M.Yamauchi h,H.Yanai ac,P.Yeh z,Y.Yuan j,Y.Yusa ar,J.Zhang ay,Z.P.Zhang ak,Y.Zheng g,and D.ˇZontar aya Aomori University,Aomori,Japanb Budker Institute of Nuclear Physics,Novosibirsk,Russiac Chiba University,Chiba,Japand Chuo University,Tokyo,Japane University of Cincinnati,Cincinnati,OH,USAf University of Frankfurt,Frankfurt,Germanyg University of Hawaii,Honolulu,HI,USAh High Energy Accelerator Research Organization(KEK),Tsukuba,Japani Hiroshima Institute of Technology,Hiroshima,Japanj Institute of High Energy Physics,Chinese Academy of Sciences,Beijing,PRChinak Institute of High Energy Physics,Vienna,Austria ℓInstitute for Theoretical and Experimental Physics,Moscow,Russiam J.Stefan Institute,Ljubljana,Slovenian Kanagawa University,Yokohama,Japano Korea University,Seoul,South Koreap Kyoto University,Kyoto,Japanq Kyungpook National University,Taegu,South Korear Institut de Physique des Hautes´Energies,Universit´e de Lausanne,Lausanne,Switzerlands University of Ljubljana,Ljubljana,Sloveniat University of Maribor,Maribor,Sloveniau University of Melbourne,Victoria,Australiav Nagoya University,Nagoya,Japanw Nara Women’s University,Nara,Japanx National Kaohsiung Normal University,Kaohsiung,Taiwany National Lien-Ho Institute of Technology,Miao Li,Taiwanz National Taiwan University,Taipei,Taiwanaa H.Niewodniczanski Institute of Nuclear Physics,Krakow,Polandab Nihon Dental College,Niigata,Japanac Niigata University,Niigata,Japanad Osaka City University,Osaka,Japanae Osaka University,Osaka,Japanaf Panjab University,Chandigarh,Indiaag Peking University,Beijing,PR Chinaah Princeton University,Princeton,NJ,USAai RIKEN BNL Research Center,Brookhaven,NY,USAaj Saga University,Saga,Japanak University of Science and Technology of China,Hefei,PR ChinaaℓSeoul National University,Seoul,South Koreaam Sungkyunkwan University,Suwon,South Koreaan University of Sydney,Sydney,NSW,Australiaao Tata Institute of Fundamental Research,Bombay,Indiaap Toho University,Funabashi,Japanaq Tohoku Gakuin University,Tagajo,Japanar Tohoku University,Sendai,Japanas University of Tokyo,Tokyo,Japanat Tokyo Institute of Technology,Tokyo,Japanau Tokyo Metropolitan University,Tokyo,Japanav Tokyo University of Agriculture and Technology,Tokyo,Japanaw Toyama National College of Maritime Technology,Toyama,Japanay University of Tsukuba,Tsukuba,Japanaz Utkal University,Bhubaneswer,Indiaba Virginia Polytechnic Institute and State University,Blacksburg,VA,USAbb Yokkaichi University,Yokkaichi,Japanbc Yonsei University,Seoul,South KoreaB events collected with the Belle detector at KEKB.Thebranching fractions B(B+→ρ0π+)=(8.0+2.3+0.7−2.0−0.7)×10−6and B(B0→ρ±π∓)=(20.8+6.0+2.8−6.3−3.1)×10−6are obtained.In addition,a90%confidence level upper limitof B(B0→ρ0π0)<5.3×10−6is reported.Key words:ρπ,branching fractionPACS:13.25.hw,14.40.Nd1on leave from Nova Gorica Polytechnic,Nova Gorica,Sloveniamodes are examined.Here and throughout the text,inclusion of charge con-jugate modes is implied and for the neutral decay,B0→ρ±π∓,the notation implies a sum over both the modes.The data sample used in this analysis was taken by the Belle detector[9]at KEKB[10],an asymmetric storage ring that collides8GeV electrons against3.5GeV positrons.This produces Υ(4S)mesons that decay into B0B pairs.The Belle detector is a general purpose spectrometer based on a1.5T su-perconducting solenoid magnet.Charged particle tracking is achieved with a three-layer double-sided silicon vertex detector(SVD)surrounded by a central drift chamber(CDC)that consists of50layers segmented into6axial and5 stereo super-layers.The CDC covers the polar angle range between17◦and 150◦in the laboratory frame,which corresponds to92%of the full centre of mass(CM)frame solid angle.Together with the SVD,a transverse momen-tum resolution of(σp t/p t)2=(0.0019p t)2+(0.0030)2is achieved,where p t is in GeV/c.Charged hadron identification is provided by a combination of three devices: a system of1188aerogelˇCerenkov counters(ACC)covering the momentum range1–3.5GeV/c,a time-of-flight scintillation counter system(TOF)for track momenta below1.5GeV/c,and dE/dx information from the CDC for particles with very low or high rmation from these three devices is combined to give the likelihood of a particle being a kaon,L K,or pion, Lπ.Kaon-pion separation is then accomplished based on the likelihood ratio Lπ/(Lπ+L K).Particles with a likelihood ratio greater than0.6are identified as pions.The pion identification efficiencies are measured using a high momentum D∗+data sample,where D∗+→D0π+and D0→K−π+.With this pion selection criterion,the typical efficiency for identifying pions in the momentum region0.5GeV/c<p<4GeV/c is(88.5±0.1)%.By comparing the D∗+data sample with a Monte Carlo(MC)sample,the systematic error in the particle identification(PID)is estimated to be1.4%for the mode with three charged tracks and0.9%for the modes with two.Surrounding the charged PID devices is an electromagnetic calorimeter(ECL) consisting of8736CsI(Tl)crystals with a typical cross-section of5.5×5.5cm2 at the front surface and16.2X0in depth.The ECL provides a photon energy resolution of(σE/E)2=0.0132+(0.0007/E)2+(0.008/E1/4)2,where E is in GeV.Electron identification is achieved by using a combination of dE/dx measure-ments in the CDC,the response of the ACC and the position and shape of the electromagnetic shower from the ECL.Further information is obtained from the ratio of the total energy registered in the calorimeter to the particle momentum,E/p lab.Charged tracks are required to come from the interaction point and have transverse momenta above100MeV/c.Tracks consistent with being an elec-tron are rejected and the remaining tracks must satisfy the pion identification requirement.The performance of the charged track reconstruction is studied using high momentumη→γγandη→π+π−π0decays.Based on the relative yields between data and MC,we assign a systematic error of2%to the single track reconstruction efficiency.Neutral pion candidates are detected with the ECL via their decayπ0→γγ. Theπ0mass resolution,which is asymmetric and varies slowly with theπ0 energy,averages toσ=4.9MeV/c2.The neutral pion candidates are selected fromγγpairs by requiring that their invariant mass to be within3σof the nominalπ0mass.To reduce combinatorial background,a selection criteria is applied to the pho-ton energies and theπ0momenta.Photons in the barrel region are required to have energies over50MeV,while a100MeV requirement is made for photons in the end-cap region.Theπ0candidates are required to have a momentum greater than200MeV/c in the laboratory frame.Forπ0s from BE2beam−p2B and the energy difference∆E=E B−E beam.Here, p B and E B are the momentum and energy of a B candidate in the CM frame and E beam is the CM beam energy.An incorrect mass hypothesis for a pion or kaon produces a shift of about46MeV in∆E,providing extra discrimination between these particles.The width of the M bc distributions is primarily due to the beam energy spread and is well modelled with a Gaussian of width 3.3MeV/c2for the modes with a neutral pion and2.7MeV/c2for the mode without.The∆E distribution is found to be asymmetric with a small tail on the lower side for the modes with aπ0.This is due toγinteractions withmaterial in front of the calorimeter and shower leakage out of the calorimeter. The∆E distribution can be well modelled with a Gaussian when no neutral particles are present.Events with5.2GeV/c2<M bc<5.3GeV/c2and|∆E|< 0.3GeV are selected for thefinal analysis.The dominant background comes from continuum e+e−→qB events and jet-like qi,j|p i||p j|P l(cosθij)i,k|p i||p k|,r l=),where L s and L qqD0π+ decays.By comparing the yields in data and MC after the likelihood ratiorequirement,the systematic errors are determined to be4%for the modes with aπ0and6%for the mode without.Thefinal variable used for continuum suppression is theρhelicity angle,θh, defined as the angle between the direction of the decay pion from theρin the ρrest frame and theρin the B rest frame.The requirement of|cosθh|>0.3 is made independently of the likelihood ratio as it is effective in suppressing the background from B decays as well as the qB events is used[14].The largest component of this background is found to come from decays of the type B→Dπ;when the D meson decays via D→π+π−,events can directly reach the signal region while the decay D→K−π+can reach the signal region with the kaon misidentified as a pion.Decays with J/ψandψ(2S) mesons can also populate the signal region if both the daughter leptons are misidentified as pions.These events are excluded by making requirements on the invariant mass of the intermediate particles:|M(π+π−)−M D0|>0.14 GeV/c2,|M(π+π0)−M D+|>0.05GeV/c2,|M(π+π−)−M J/ψ|>0.07GeV/c2 and|M(π+π−)−Mψ(2S)|>0.05GeV/c2.The widest cut is made around the D0mass to account for the mass shift due to misidentifying the kaons in D0 decays as pions.Fig.1shows the∆E and M bc distributions for the three modes analysed after all the selection criteria have been applied.The∆E and M bc plots are shown for events that lie within3σof the nominal M bc and∆E values,respectively. The signal yields are obtained by performing maximum likelihoodfits,each using a single signal function and one or more background functions.The signal functions are obtained from the MC and adjusted based on comparisons of B+→B0are assumed to be equal.The M bc distribution for all modes isfitted with a single Gaussian and an ARGUS background function[15].The normalization of the ARGUS function is left tofloat and shape of the function isfixed from the∆E sideband:−0.25 GeV<∆E<−0.08GeV and5.2GeV/c2<M bc<5.3GeV/c2.For the mode with only charged pions in thefinal state,the∆E distribution isfitted with a single Gaussian for the signal and a linear function withfixed shape for the continuum background.The normalization of the linear function is left to float and the slope isfixed from the M bc sideband,5.2GeV/c2<M bc<5.26GeV/c2,|∆E|<0.3GeV.There are also other rare B decays that are expected to contaminate the∆E distribution.For the mode without aπ0,these modes are of the type B0→h+h−(where h denotes aπor K),B→ρρ(including all combinations of charged and neutralρmesons,where the polarizations of theρmesons are assumed to be longitudinal)and B→Kππ(including the decays B+→ρ0K+,B+→K∗0π+,B+→K∗0(1430)0π+,B+→f0(980)K+ and B+→f0(1370)K+)[16].These background modes are accounted for by using smoothed histograms whose shapes have been determined by combining MC distributions.The three B→ρρmodes are combined into one histogram. The normalization of this component is allowed tofloat in thefit due to the uncertainty in the branching fractions of the B→ρρmodes.Likewise,the B→hh and all the B→Kππmodes are combined to form one hh and one Kππcomponent.The normalizations of these components arefixed to their expected yields,which are calculated using efficiencies determined by MC and branching fractions measured by previous Belle analyses[16,17].The∆Efits for the modes with aπ0in thefinal state have the signal compo-nent modelled by a Crystal Ball function[18]to account for the asymmetry in the∆E distribution.As for the B+→ρ0π+mode,the continuum background is modelled by a linear function withfixed slope.Unlike the B+→ρ0π+mode, a component is included for the background from the b→c transition.The pa-rameterization for rare B decays includes one component for the B→Kππ0 modes(B0→ρ+K−and B0→K∗+π−)[19]and one for all the B→ρρmodes.The normalization of the B→ρρcomponent is left tofloat while the other components from B decays arefixed to their expected yields.Table1summarizes the results of the∆Efits,showing the number of events, signal yields,reconstruction efficiencies,statistical significance and branching fractions or upper limits for eachfit.The statistical significance is defined assystematic error in thefitted signal yield is estimated by independently varying eachfixed parameter in thefit by1σ.Thefinal results are B(B+→ρ0π+)=(8.0+2.3+0.7−2.0−0.7)×10−6and B(B0→ρ±π∓)=(20.8+6.0+2.8−6.3−3.1)×10−6where thefirst error is statistical and the second is systematic.For theρ0π0mode,one standard deviation of the systematic error is added to the statistical limit to obtain a conservative upper limit at90%confidence of5.3×10−6.The possibility of a nonresonant B→πππbackground is also examined.To check for this type of background,the M bc and∆E yields are determined for differentππinvariant mass bins.Byfitting the M bc distribution inππinvariant mass bins with B→ρπand B→πππMC distributions,the nonresonant contribution is found to be below4%.To account for this possible background, errors3.7%and3.2%are added in quadrature to the systematic errors of the ρ+π−andρ0π+modes,respectively.Theππinvariant mass distributions are shown in Fig.2.Two plots are shown for theρ+π−andρ0π+modes,one with events from the M bc sideband superimposed over the events from the signal region(upper)and one with events from signal MC superimposed over events from the signal region with the sideband subtracted(lower).Fig.3 shows the distribution of the helicity variable,cosθh,for the two modes with all selection criteria applied except the helicity condition.Events fromρπdecays are expected to follow a cos2θdistribution while nonresonant and other background decays have an approximately uniform distribution.The helicity plots are obtained byfitting the M bc distribution in eight helicity bins ranging from−1to1.The M bc yield is then plotted against the helicity bin for each mode and the expected MC signal distributions are superimposed.Both the ππmass spectrum and the helicity distributions provide evidence that the signal events are consistent with being fromρπdecays.The results obtained here can be used to calculate the ratio of branching frac-tions R=B(B0→ρ±π∓)/B(B+→ρ0π+),which gives R=2.6±1.0±0.4, where thefirst error is statistical and second is systematic.This is consistent with values obtained by CLEO[20]and BaBar[21,22]as shown in Table2. Theoretical calculations done at tree level assuming the factorization approx-imation for the hadronic matrix elements give R∼6[3].Calculations that include penguin contributions,off-shell B∗excited states or additionalππres-onances[4–8]might yield better agreement with the the measured value of R.In conclusion,statistically significant signals have been observed in the B→ρπmodes using a31.9×106BWe wish to thank the KEKB accelerator group for the excellent operation of the KEKB accelerator.We acknowledge support from the Ministry of Ed-ucation,Culture,Sports,Science,and Technology of Japan and the Japan Society for the Promotion of Science;the Australian Research Council and the Australian Department of Industry,Science and Resources;the National Science Foundation of China under contract No.10175071;the Department of Science and Technology of India;the BK21program of the Ministry of Education of Korea and the CHEP SRC program of the Korea Science and Engineering Foundation;the Polish State Committee for Scientific Research under contract No.2P03B17017;the Ministry of Science and Technology of the Russian Federation;the Ministry of Education,Science and Sport of the Republic of Slovenia;the National Science Council and the Ministry of Education of Taiwan;and the U.S.Department of Energy.References[1] A.E.Snyder and H.R.Quinn,Phys.Rev.D48,2139(1993).[2]I.Bediaga,R.E.Blanco,C.G¨o bel,and R.M´e ndez-Galain,Phys.Rev.Lett.81,4067(1998).[3]M.Bauer,B.Stech,and M.Wirbel,Z.Phys.C34,103(1987).[4] A.Deandrea et al.,Phys.Rev.D62,036001(2000).[5]Y.H.Chen,H.Y.Cheng,B.Tseng and K.C.Yang,Phys.Rev.D60,094014(1999).[6] C.D.Lu and M.Z.Yang,Eur.Phys.J C23,275(2002).[7]J.Tandean and S.Gardner,SLAC-PUB-9199;hep-ph/0204147.[8]S.Gardner and Ulf-G.Meißner,Phys.Rev.D65,094004(2002).[9]Belle Collaboration,A.Abashian et al.,Nucl.Instr.and Meth.A479,117(2002).[10]E.Kikutani ed.,KEK Preprint2001-157(2001),to appear in Nucl.Instr.andMeth.A.[11]G.C.Fox and S.Wolfram,Phys.Rev.Lett.41,1581(1978).[12]This modification of the Fox-Wolfram moments wasfirst proposed in a seriesof lectures on continuum suppression at KEK by Dr.R.Enomoto in May and June of1999.For a more detailed description see Belle Collaboration,K.Abe et al.,Phys.Lett.B511,151(2001).[13]CLEO Collaboration,D.M.Asner et al.,Phys.Rev.D53,1039(1996).[14]These MC events are generated with the CLEO group’s QQ program,see/public/CLEO/soft/QQ.The detector response is simulated using GEANT,R.Brun et al.,GEANT 3.21,CERN Report DD/EE/84-1,1984.[15]The ARGUS Collaboration,H.Albrecht et al.,Phys.Lett.B241,278(1990).[16]Belle Collaboration,A.Garmash et al.,Phys.Rev.D65,092005(2002).[17]Belle Collaboration,K.Abe et al.,Phys.Rev.Lett.87,101801(2001).[18]J.E.Gaiser et al.,Phys.Rev.D34,711(1986).[19]Belle Collaboration,K.Abe et al.,BELLE-CONF-0115,submitted as acontribution paper to the2001International Europhysics Conference on High Energy Physics(EPS-HEP2001).[20]CLEO Collaboration,C.P.Jessop et al.,Phys.Rev.Lett.85,2881(2000).[21]Babar Collaboration,B.Aubert et al.,submitted as a contribution paper tothe20th International Symposium on Lepton and Photon Interactions at High Energy(LP01);hep-ex/0107058.[22]BaBar Collaboration,B.Aubert et al.,submitted as a contribution paper tothe XXXth International Conference on High Energy Physics(ICHEP2000);hep-ex/0008058.Table1∆Efit results.Shown for each mode are the number of events in thefit,the signal yield,the reconstruction efficiency,the branching fraction(B)or90%confidence level upper limit(UL)and the statistical significance of thefit.Thefirst error in the branching fraction is statistical,the second is systematic.ρ0π+15424.3+6.9−6.29.68.0+2.3+0.7−2.0−0.74.4σρ+π−30144.6+12.8−13.46.820.8+6.0+2.8−6.3−3.13.7σρ0π0116−4.4±8.58.5<5.3-Experiment B(B0→ρ±π∓)×10−6B(B+→ρ0π+)×10−6RE v e n t s /16 M e VE v e n t s /3 M e V /c2(b) ρ0π+Signal backgrd02.557.51012.51517.52022.55.25.225 5.25 5.2755.3E v e n t s /18 M e VE v e n t s /2 M e V /c2(d) ρ+π-Signal backgrd051015202530355.25.225 5.25 5.2755.3∆E(GeV)E v e n t s /18 M e V(e) ρ0π024681012-0.2-0.10.10.2(GeV/c 2)E v e n t s /2 M e V /c2M bc (f) ρ0πSignal backgrd02468101214165.25.225 5.25 5.2755.3Fig.1.The ∆E (left)and M bc (right)fits to the three B →ρπmodes:ρ0π+,ρ+π−and ρ0π0.The histograms show the data,the solid lines show the total fit and the dashed lines show the continuum component.In (a)the contribution from the B →ρρand B →hh modes is shown by the cross hatched component.In (c)the cross hatched component shows the contribution from the b →c transition and B →ρρmodes.102030405060+0(GeV/c 2)E v e n t s /0.1 G e V /c2M(π+π0)(GeV/c 2)E v e n t s /0.1 G e V /c2(GeV/c 2)E v e n t s /0.1 G e V /c2+-(GeV/c 2)E v e n t s /0.1 G e V /c2M(π+ π-)510152025Fig.2.The M (ππ)distributions for B 0→ρ±π∓(left)and B +→ρ0π+(right)events in the signal region.Plots (a)and (b)show sideband events superimposed;plots (c)and (d)show the sideband subtracted plots with signal MC superimposed.-1-0.500.51M b c y i e l d (E v e n t s )cos θh-1-0.500.51M b c y i e l d (E v e n t s )cos θhFig.3.The ρmeson helicity distributions for B 0→ρ±π∓(a)and B +→ρ0π+(b).Signal MC distributions are shown superimposed.。

FINDINGS FROM ROUTINE LABORATORY TESTS
• • Leukocyte count and sedimentation rate may be normal or elevated • • Thrombocytopenia • • Hyponatremia • • Elevated serum transaminases • • Elevated creatine phosphokinase • • Myoglobinuria • • Sterile pyuria
commonly used serologic tests lack precision and accuracy
cross-reactivity with syphilis, relapsing fever, leptopsirosis, HIV, RMSF, infectious mononucleosis,
SIGNS/SYMPTOMS
LOCALIZED (WITHIN 3-30 DAYS) 1. Erythema migrans (EM) — red ring-like or homogenous expanding rash 2. May have flu-like symptoms
DISSEMINATED (WITHIN WEEKS TO MONTHS) • Fatigue, Malaise ,Lethargy, Headache ,Stiff neck,Myalgia, Arthralgia ,Regional or generalized lymphadenopathy, Fever,Erythema migrans (single or multiple lesion)
DIAGNOSTIC LABORATORY CRITERIA

Chapter7Ordinary Differential EquationsMatlab has several different functions for the numerical solution of ordinary dif-ferential equations.This chapter describes the simplest of these functions and then compares all of the functions for efficiency,accuracy,and special features.Stiffness is a subtle concept that plays an important role in these comparisons.7.1Integrating Differential EquationsThe initial value problem for an ordinary differential equation involvesfinding a function y(t)that satisfiesdy(t)=f(t,y(t))dttogether with the initial conditiony(t0)=y0A numerical solution to this problem generates a sequence of values for the indepen-dent variable,t0,t1,...,and a corresponding sequence of values for the dependent variable,y0,y1,...,so that each y n approximates the solution at t n y n≈y(t n),n=0,1,...Modern numerical methods automatically determine the step sizesh n=t n+1−t nso that the estimated error in the numerical solution is controlled by a specified tolerance.The Fundamental Theorem of Calculus gives us an important connection be-tween differential equations and integrals.t+hy(t+h)=y(t)+f(s,y(s))dst12Chapter7.Ordinary Differential Equations We cannot use numerical quadrature directly to approximate the integral because we do not know the function y(s)and so cannot evaluate the integrand.Nevertheless, the basic idea is to choose a sequence of values of h so that this formula allows us to generate our numerical solution.One special case to keep in mind is the situation where f(t,y)is a function of t alone.The numerical solution of such simple differential equations is then just a sequence of quadratures.y n+1=y n+ tn+1t nf(s)dsThrough this chapter,we frequently use“dot”notation for derivatives,˙y=dy(t)dtand¨y=d2y(t)dt27.2Systems of EquationsMany mathematical models involve more than one unknown function,and second and higher order derivatives.These models can be handled by making y(t)a vector-valued function of t.Each component is either one of the unknown functions or one of its derivatives.The Matlab vector notation is particularly convenient here.For example,the second-order differential equation describing a simple har-monic oscillator¨x(t)=−x(t)becomes twofirst-order equations.The vector y(t)has two components,x(t)and itsfirst derivative˙x(t),y(t)=x(t)˙x(t)Using this vector,the differential equation is˙y(t)=˙x(t)−x(t)=y2(t)−y1(t)The Matlab function defining the differential equation has t and y as input arguments and should return f(t,y)as a column vector.For the harmonic oscillator, the function isfunction ydot=harmonic(t,y)ydot=[y(2);-y(1)]A fancier version uses matrix multiplication in an inline function.f=inline(’[01;-10]*y’,’t’,’y’);7.3.Linearized Differential Equations3 In both cases,the variable t has to be included as thefirst argument,even though it is not explicitly involved in the differential equation.A slightly more complicated example,the two-body problem,describes the orbit of one body under the gravitational attraction of a much heavier ing Cartesian coordinates,u(t)and v(t),centered in the heavy body,the equations are ¨u(t)=−u(t)/r(t)3¨v(t)=−v(t)/r(t)3wherer(t)=u(t)+v(t)The vector y(t)has four components,y(t)=u(t)v(t)˙u(t)˙v(t)The differential equation is˙y(t)=˙u(t)˙v(t)−u(t)/r(t)3−v(t)/r(t)3The Matlab function could befunction ydot=twobody(t,y)r=sqrt(y(1)^2+y(2)^2);ydot=[y(3);y(4);-y(1)/r^3;-y(2)/r^3];A more compact Matlab function isfunction ydot=twobody(t,y)ydot=[y(3:4);-y(1:2)/norm(y(1:2))^3]Despite the use of vector operations,the second M-file is not significantly more efficient than thefirst.7.3Linearized Differential EquationsThe local behavior of the solution to a differential equation near any point(t c,y c) can be analyzed by expanding f(t,y)in a two-dimensional Taylor series.f(t,y)=f(t c,y c)+α(t−t c)+J(y−y c)+...whereα=∂f∂t(t c,y c),J=∂f∂y(t c,y c)4Chapter7.Ordinary Differential Equations The most important term in this series is usually the one involving J,the Jacobian. For a system of differential equations with n components,d dty1(t)y2(t)...y n(t)=f1(t,y1,...,y n)f2(t,y1,...,y n)...f n(t,y1,...,y n)the Jacobian is an n-by-n matrix of partial derivativesJ=∂f1∂y1∂f1∂y2...∂f1∂y n ∂f2∂y1∂f2∂y2...∂f2∂y n.........∂f n∂y1∂f n∂y2...∂f n∂y nThe influence of the Jacobian on the local behavior is determined by the solution to the linear system of ordinary differential equations˙y=JyLetλk=µk+iνk be the eigenvalues of J andΛ=diag(λk)the diagonal eigenvalue matrix.If there is a linearly independent set of corresponding eigenvectors V,then J=VΛV−1The linear transformationV x=ytransforms the local system of equations into a set of decoupled equations for the individual components of x,˙x k=λk x kThe solutions arex k(t)=eλk(t−t c)x(t c)A single component x k(t)grows with t ifµk is positive,decays ifµk is negative, and oscillates ifνk is nonzero.The components of the local solution y(t)are linear combinations of these behaviors.For example the harmonic oscillator˙y=01−10yis a linear system.The Jacobian is simply the matrixJ=01−107.4.Single-Step Methods 5The eigenvalues of J are ±i and the solutions are purely oscillatory linear combi-nations of e it and e −it .A nonlinear example is the two-body problem ˙y (t )= y 3(t )y 4(t )−y 1(t )/r (t )3−y 2(t )/r (t )3wherer (t )= y 1(t )2+y 2(t )2In an exercise,we ask you to show that the Jacobian for this system is J =1r 5 00r 50000r 52y 21−y 223y 1y 2003y 1y 22y 22−y 2100It turns out that the eigenvalues of J just depend on the radius r (t )λ=1r3/2 √2i −√2−iWe see that one eigenvalue is real and positive,so the corresponding component of the solution is growing.One eigenvalue is real and negative,corresponding to a decaying component.Two eigenvalues are purely imaginary,corresponding to os-cillatory components.However,the overall global behavior of this nonlinear system is quite complicated,and is not described by this local linearized analysis.7.4Single-Step MethodsThe simplest numerical method for the solution of initial value problems is Euler’s method.It uses a fixed step size h and generates the approximate solution byy n +1=y n +hf (t n ,y n )t n +1=t n +hThe Matlab code would use an initial point t0,a final point tfinal ,an initial value y0,a step size h ,and an inline function or function handle f .The primary loop would simply bet =t0;y =y0;while t <=tfinaly =y +h*feval(f,t,y)t =t +hend6Chapter7.Ordinary Differential Equations Note that this works perfectly well if y0is a vector and f returns a vector.As a quadrature rule for integrating f(t),Euler’s method corresponds to a rectangle rule where the integrand is evaluated only once,at the left hand endpoint of the interval.It is exact if f(t)is constant,but not if f(t)is linear.So the error is proportional to h.Tiny steps are needed to get even a few digits of accuracy. But,from our point of view,the biggest defect of Euler’s method is that it does not provide an error estimate.There is no automatic way to determine what step size is needed to achieve a specified accuracy.If Euler’s method is followed by a second function evaluation,we begin to get a viable algorithm.There are two natural possibilities,corresponding to the midpoint rule and the trapezoid rule for quadrature.The midpoint analog uses Euler to step halfway across the interval,evaluates the function at this intermediate point,then uses that slope to take the actual step.s1=f(t n,y n)s2=f(t n+h2,y n+h2s1)y n+1=y n+hs2t n+1=t n+hThe trapezoid analog uses Euler to take a tentative step across the interval,evaluates the function at this exploratory point,then averages the two slopes to take the actual step.s1=f(t n,y n)s2=f(t n+h,y n+hs1)y n+1=y n+h s1+s22t n+1=t n+hIf we were to use both of these methods simultaneously,they would produce two different values for y n+1.The difference between the two values would provide an error estimate and a basis for picking the step size.Furthermore,an extrapolated combination of the two values would be more accurate than either one individually.Continuing with this approach is the idea behind single-step methods for inte-grating ODEs.The function f(t,y)is evaluated several times for values of t between t n and t n+1and values of y obtained by adding linear combinations of the values of f to y n.The actual step is taken using another linear combination of the function values.Modern versions of single-step methods use yet another linear combination of function values to estimate error and determine step size.Single-step methods are often called Runge-Kutta methods,after the two Ger-man applied mathematicians whofirst wrote about them around1905.The classical Runge-Kutta method was widely used for hand computation before the invention of digital computers and is still popular today.It uses four function evaluations per step.s1=f(t n,y n)7.4.Single-Step Methods7s2=f(t n+h2,y n+h2s1)s3=f(t n+h2,y n+h2s2)s4=f(t n+h,y n+hs3)y n+1=y n+h6(s1+2s2+2s3+s4)t n+1=t n+hIf f(t,y)does not depend on y,then classical Runge-Kutta has s2=s3and the method reduces to Simpson’s quadrature rule.Classical Runge-Kutta does not provide an error estimate.The method is sometimes used with a step size h and again with step size h/2to obtain an error estimate,but we now know more efficient methods.Several of the ODE solvers in Matlab,including the textbook solver we describe later in this chapter,are single-step or Runge-Kutta solvers.A general single-step method is characterized by a number of parameters,αi,βi,j,γi,andδi. There are k stages.Each stage computes a slope,s i,by evaluating f(t,y)for a particular value of t and a value of y obtained by taking linear combinations of the previous slopes.s i=f(t n+αi h,y n+h i−1j=1βi,j s j),i=1,...,kThe proposed step is also a linear combination of the slopes.y n+1=y n+hki=1γi s iAn estimate of the error that would occur with this step is provided by yet another linear combination of the slopes.e n+1=hki=1δi s iIf this error is less than the specified tolerance,then the step is successful and y n+1 is accepted.If not,the step is a failure and y n+1is rejected.In either case,the error estimate is used to compute the step size h for the next step.The parameters in these methods are determined by matching terms in Taylor series expansions of the slopes.These series involve powers of h and products of various partial derivatives of f(t,y).The order of a method is the exponent of the smallest power of h that cannot be matched.It turns out that one,two,three,or four stages yield methods of order one,two,three or four,respectively.But it takes six stages to obtain afifth-order method.The classical Runge-Kutta method has four stages and is fourth-order.The names of the Matlab ODE solvers are all of the form odennxx with digits nn indicating the order of the underlying method and a possibly empty xx indicating8Chapter7.Ordinary Differential Equationssome special characteristic of the method.If the error estimate is obtained by comparing formulas with different orders,the digits nn indicate these orders.For example,ode45obtains its error estimate by comparing a fourth-order and afifth-order formula.7.5The BS23AlgorithmOur textbook function ode23tx is a simplified version of the function ode23that is included with Matlab.The algorithm is due to Bogacki and Shampine[3,6].The “23”in the function names indicates that two simultaneous single step formulas, one of second order and one of third order,are involved.The method has three stages,but there are four slopes s i because,after the first step,the s1for one step is the s4from the previous step.The essentials are s1=f(t n,y n)s2=f(t n+h2,y n+h2s1)s3=f(t n+34h,y n+34hs2)t n+1=t n+hy n+1=y n+h9(2s1+3s2+4s3)s4=f(t n+1,y n+1)e n+1=h72(−5s1+6s2+8s3−9s4)The simplified pictures infigure7.1show the starting situation and the three stages.We start at a point(t n,y n)with an initial slope s1=f(t n,y n)and an estimate of a good step size,h.Our goal is to compute an approximate solution y n+1at t n+1=t n+h that agrees with true solution y(t n+1)to within the specified tolerances.Thefirst stage uses the initial slope s1to take an Euler step halfway across the interval.The function is evaluated there to get the second slope,s2.This slope is used to take an Euler step three quarters of the way across the interval.The function is evaluated again to get the third slope,s3.A weighted average of the three slopess=19(2s1+3s2+4s3)is used for thefinal step all the way across the interval to get a tentative value for y n+1.The function is evaluated once more to get s4.The error estimate then uses all four slopes.e n+1=h72(−5s1+6s2+8s3−9s4)If the error is within the specified tolerance,then the step is successful,the tentative value of y n+1is accepted,and s4becomes the s1of the next step.If the error isFigure7.1.BS23algorithmtoo large,then the tentative y n+1is rejected and the step must be redone.In either case,the error estimate e n+1provides the basis for determining the step size h for the next step.Thefirst input argument of ode23tx specifies the function f(t,y).This argu-ment can take any of three different forms:•A character string involving t and/or y•An inline function•A function handleThe inline function or the M-file should accept two arguments,usually,but not necessarily,t and y.The result of evaluating the character string or the function should be a column vector containing the values of the derivatives,dy/dt.The second input argument of ode23tx is a vector,tspan,with two compo-nents,t0and tfinal.The integration is carried out over the interval t0≤t≤t finalOne of the simplifications in our textbook code is this form of tspan.Other Matlab ODE solvers allow moreflexible specifications of the integration interval.The third input argument is a column vector,y0,providing the initial value of y0=y(t0).The length of y0tells ode23tx the number of differential equations in the system.10Chapter7.Ordinary Differential EquationsA fourth input argument is optional,and can take two different forms.The simplest,and most common,form is a scalar numerical value,rtol,to be used as the relative error tolerance.The default value for rtol is10−3,but you can provide a different value if you want more or less accuracy.The more complicated possibility for this optional argument is the structure generated by the Matlab function odeset.This function takes pairs of arguments that specify many different options for the Matlab ODE solvers.For ode23tx,you can change the default values of three quantities,the relative error tolerance,the absolute error tolerance, and the M-file that is called after each successful step.The statement opts=odeset(’reltol’,1.e-5,’abstol’,1.e-8,...’outputfcn’,@myodeplot)creates a structure that specifies the relative error tolerance to be10−5,the absolute error tolerance to be10−8,and the output function to be myodeplot.The output produced by ode23tx can be either graphic or numeric.With no output arguments,the statementode23tx(F,tspan,y0);produces a dynamic plot of all the components of the solution.With two output arguments,the statement[tout,yout]=ode23tx(F,tspan,y0);generates a table of values of the solution.7.6ode23txLet’s examine the code for ode23tx.Here is the preamble.function[tout,yout]=ode23tx(F,tspan,y0,arg4,varargin)%ODE23TX Solve non-stiff differential equations.%Textbook version of ODE23.%%ODE23TX(F,TSPAN,Y0)with TSPAN=[T0TFINAL]integrates%the system of differential equations y’=f(t,y)from%t=T0to t=TFINAL.The initial condition is%y(T0)=Y0.The input argument F is the name of an%M-file,or an inline function,or simply a character%string,defining f(t,y).This function must have two%input arguments,t and y,and must return a column vector%of the derivatives,y’.%%With two output arguments,[T,Y]=ODE23TX(...)returns%a column vector T and an array Y where Y(:,k)is the%solution at T(k).%7.6.ode23tx11%With no output arguments,ODE23TX plots the solution.%%ODE23TX(F,TSPAN,Y0,RTOL)uses the relative error%tolerance RTOL instead of the default 1.e-3.%%ODE23TX(F,TSPAN,Y0,OPTS)where OPTS=ODESET(’reltol’,...%RTOL,’abstol’,ATOL,’outputfcn’,@PLTFN)uses relative error %RTOL instead of 1.e-3,absolute error ATOL instead of% 1.e-6,and calls PLTFN instead of ODEPLOT after each step.%%More than four input arguments,ODE23TX(F,TSPAN,Y0,RTOL,%P1,P2,..),are passed on to F,F(T,Y,P1,P2,..).%%ODE23TX uses the Runge-Kutta(2,3)method of Bogacki and%Shampine.%%Example%tspan=[02*pi];%y0=[10]’;%F=’[01;-10]*y’;%ode23tx(F,tspan,y0);%%See also ODE23.Here is the code that parses the arguments and initializes the internal variables.rtol= 1.e-3;atol= 1.e-6;plotfun=@odeplot;if nargin>=4&isnumeric(arg4)rtol=arg4;elseif nargin>=4&isstruct(arg4)if~isempty(arg4.RelTol),rtol=arg4.RelTol;endif~isempty(arg4.AbsTol),atol=arg4.AbsTol;endif~isempty(arg4.OutputFcn),plotfun=arg4.OutputFcn;end endt0=tspan(1);tfinal=tspan(2);tdir=sign(tfinal-t0);plotit=(nargout==0);threshold=atol/rtol;hmax=abs(0.1*(tfinal-t0));t=t0;y=y0(:);%Make F callable by feval.12Chapter7.Ordinary Differential Equations if ischar(F)&exist(F)~=2F=inline(F,’t’,’y’);elseif isa(F,’sym’)F=inline(char(F),’t’,’y’);end%Initialize output.if plotitfeval(plotfun,tspan,y,’init’);elsetout=t;yout=y.’;endThe computation of the initial step size is a delicate matter because it requires some knowledge of the overall scale of the problem.s1=feval(F,t,y,varargin{:});r=norm(s1./max(abs(y),threshold),inf)+realmin;h=tdir*0.8*rtol^(1/3)/r;Here is the beginning of the main loop.The integration starts at t=t0and increments t until it reaches t final.It is possible to go“backward,”that is,have t final<t0.while t~=tfinalhmin=16*eps*abs(t);if abs(h)>hmax,h=tdir*hmax;endif abs(h)<hmin,h=tdir*hmin;end%Stretch the step if t is close to tfinal.if 1.1*abs(h)>=abs(tfinal-t)h=tfinal-t;endHere is the actual computation.Thefirst slope s1has already been computed.The function defining the differential equation is evaluated three more times to obtain three more slopes.s2=feval(F,t+h/2,y+h/2*s1,varargin{:});s3=feval(F,t+3*h/4,y+3*h/4*s2,varargin{:});tnew=t+h;ynew=y+h*(2*s1+3*s2+4*s3)/9;s4=feval(F,tnew,ynew,varargin{:});7.6.ode23tx13 Here is the error estimate.The norm of the error vector is scaled by the ratio of the absolute tolerance to the relative tolerance.The use of the smallestfloating-point number,realmin,prevents err from being exactly zero.e=h*(-5*s1+6*s2+8*s3-9*s4)/72;err=norm(e./max(max(abs(y),abs(ynew)),threshold),...inf)+realmin;Here is the test to see if the step is successful.If it is,the result is plotted or appended to the output vector.If it is not,the result is simply forgotten.if err<=rtolt=tnew;y=ynew;if plotitif feval(plotfun,t,y,’’);breakendelsetout(end+1,1)=t;yout(end+1,:)=y.’;ends1=s4;%Reuse final function value to start new step.endThe error estimate is used to compute a new step size.The ratio rtol/err is greater than one if the current step is successful,or less than one if the current step fails.A cube root is involved because the BS23is a third-order method.This means that changing tolerances by a factor of eight will change the typical step size,and hence the total number of steps,by a factor of two.The factors0.8and5prevent excessive changes in step size.%Compute a new step size.h=h*min(5,0.8*(rtol/err)^(1/3));Here is the only place where a singularity would be detected.if abs(h)<=hminwarning(sprintf(...’Step size%e too small at t=%e.\n’,h,t));t=tfinal;endendThat ends the main loop.The plot function might need tofinish its work.if plotitfeval(plotfun,[],[],’done’);end14Chapter7.Ordinary Differential Equations 7.7ExamplesPlease sit down in front of a computer running Matlab.Make sure ode23tx is in your current directory or on your Matlab path.Start your session by entering ode23tx(’0’,[010],1)This should produce a plot of the solution of the initial value problem dy=0dty(0)=10≤t≤10The solution,of course,is a constant function,y(t)=1.Now you can press the up arrow key,use the left arrow key to space over to the’0’,and change it to something more interesting.Here are some examples.At first,we’ll change just the’0’and leave the[010]and1alone.F Exact solution’0’1’t’1+t^2/2’y’exp(t)’-y’exp(-t)’1/(1-3*t)’1-log(1-3*t)/3(Singular)’2*y-y^2’2/(1+exp(-2*t))Make up some of your own examples.Change the initial condition.Change the accuracy by including1.e-6as the fourth argument.Now let’s try the harmonic oscillator,a second-order differential equation writ-ten as a pair of twofirst-order equations.First,create an inline function to specify the e eitherF=inline(’[y(2);-y(1)]’,’t’,’y’)orF=inline(’[01;-10]*y’,’t’,’y’)Then,the statementode23tx(F,[02*pi],[1;0])plots two functions of t that you should recognize.If you want to produce a phase plane plot,you have two choices.One possibility is to capture the output and plot it after the computation is complete.[t,y]=ode23tx(F,[02*pi],[1;0])plot(y(:,1),y(:,2),’-o’)axis([-1.2 1.2-1.2 1.2])axis square7.7.Examples15The more interesting possibility is to use a function that plots the solution while it is being computed.Matlab provides such a function in odephas2.m.It is accessed by using odeset to create an options structureopts=odeset(’reltol’,1.e-4,’abstol’,1.e-6,...’outputfcn’,@odephas2);If you want to provide your own plotting function,it should be something like function flag=phaseplot(t,y,job)persistent pif isequal(job,’init’)p=plot(y(1),y(2),’o’,’erasemode’,’none’);axis([-1.21.2-1.2 1.2])axis squareflag=0;elseif isequal(job,’’)set(p,’xdata’,y(1),’ydata’,y(2))drawnowflag=0;endThis is withopts=odeset(’reltol’,1.e-4,’abstol’,1.e-6,...’outputfcn’,@phaseplot);Once you have decided on a plotting function and created an options structure,you can compute and simultaneously plot the solution withode23tx(F,[02*pi],[1;0],opts)Try this with other values of the tolerances.Issue the command type twobody to see if there is an M-file twobody.m on your path.If not,find the three lines of code earlier in this chapter and create your own M-file.Then tryode23tx(@twobody,[02*pi],[1;0;0;1]);The code,and the length of the initial condition,indicate that the solution has four components.But the plot shows only three.Why?Hint:find the zoom button on thefigure window toolbar and zoom in on the blue curve.You can vary the initial condition of the two-body problem by changing the fourth component.y0=[1;0;0;change_this];ode23tx(@twobody,[02*pi],y0);Graph the orbit,and the heavy body at the origin,with16Chapter7.Ordinary Differential Equations y0=[1;0;0;change_this];[t,y]=ode23tx(@twobody,[02*pi],y0);plot(y(:,1),y(:,2),’-’,0,0,’ro’)axis equalYou might also want to use something other than2πfor tfinal.7.8Lorenz AttractorOne of the world’s most extensively studied ordinary differential equations is the Lorenz chaotic attractor.It wasfirst described in1963by Edward Lorenz,an MIT mathematician and meteorologist,who was interested influidflow models of the earth’s atmosphere.An excellent reference is a book by Colin Sparrow[8].We have chosen to express the Lorenz equations in a somewhat unusual way, involving a matrix-vector product.˙y=AyThe vector y has three components that are functions of t,y(t)=y1(t)y2(t)y3(t)Despite the way we have written it,this is not a linear system of differential equa-tions.Seven of the nine elements in the3-by-3matrix A are constant,but the other two depend on y2(t).A=−β0y20−σσ−y2ρ−1Thefirst component of the solution,y1(t),is related to the convection in the atmo-sphericflow,while the other two components are related to horizontal and vertical temperature variation.The parameterσis the Prandtl number,ρis the normal-ized Rayleigh number,andβdepends on the geometry of the domain.The most popular values of the parameters,σ=10,ρ=28,andβ=8/3,are outside the ranges associated with the earth’s atmosphere.The deceptively simple nonlinearity introduced by the presence of y2in the system matrix A changes everything.There are no random aspects to these equa-tions,so the solutions y(t)are completely determined by the parameters and the initial conditions,but their behavior is very difficult to predict.For some values of the parameters,the orbit of y(t)in three-dimensional space is known as a strange attractor.It is bounded,but not periodic and not convergent.It never intersects itself.It ranges chaotically back and forth around two different points,or attractors. For other values of the parameters,the solution might converge to afixed point, diverge to infinity,or oscillate periodically.Figure 7.2.Three components of LorenzattractorFigure 7.3.Phase plane plot of Lorenz attractorLet’s think of η=y 2as a free parameter,restrict ρto be greater than one,and study the matrixA = −β0η0−σσ−ηρ−1It turns out that A is singular if and only if η=±β(ρ−1)18Chapter 7.Ordinary Differential EquationsThe corresponding null vector,normalized so that its second component is equal to η,isρ−1ηηWith two different signs for η,this defines two points in three-dimensional space.These points are fixed points for the differential equation.Ify (t 0)= ρ−1ηηthen,for all t ,˙y (t )=00and so y (t )never changes.However,these points are unstable fixed points.If y (t )does not start at one of these points,it will never reach either of them;if it tries to approach either point,it will be repulsed.We have provided an M-file,lorenzgui.m ,that facilitates experiments with the Lorenz equations.Two of the parameters,β=8/3and σ=10,are fixed.A uicontrol offers a choice among several different values of the third parameter,ρ.A simplified version of the program for ρ=28would begin withrho =28;sigma =10;beta =8/3;eta =sqrt(beta*(rho-1));A =[-beta 0eta0-sigma sigma -eta rho -1];The initial condition is taken to be near one of the attractors.yc =[rho-1;eta;eta];y0=yc +[0;0;3];The time span is infinite,so the integration will have to be stopped by another uicontrol .tspan =[0Inf];opts =odeset(’reltol’,1.e-6,’outputfcn’,@lorenzplot);ode45(@lorenzeqn,tspan,y0,opts,A);The matrix A is passed as an extra parameter to the integrator,which sends it on to lorenzeqn ,the subfunction defining the differential equation.The extra parameter machinery included in the function functions allows lorenzeqn to be written in a particularly compact manner.。

Java Lambda 表达式是Java 8 引入的一个新的功能，可以说是模拟函数式编程的一个语法糖，类似于Javascript 中的闭包，但又有些不同，主要目的是提供一个函数化的语法来简化我们的编码。

Lambda 基本语法Lambda 的基本结构为(arguments) -> body，有如下几种情况：∙参数类型可推导时，不需要指定类型，如(a) -> System.out.println(a)∙当只有一个参数且类型可推导时，不强制写(), 如a -> System.out.println(a)∙参数指定类型时，必须有括号，如(int a) -> System.out.println(a)∙参数可以为空，如() -> Sys tem.out.println(“hello”)∙body 需要用{}包含语句，当只有一条语句时{}可省略常见的写法如下：函数式接口FunctionalInterface概念Java Lambda 表达式以函数式接口为基础。

什么是函数式接口（FunctionalInterface）？简单说来就是只有一个方法（函数）的接口，这类接口的目的是为了一个单一的操作，也就相当于一个单一的函数了。

常见的接口如：Runnable, Comparator 都是函数式接口，并且都标注了注解@FunctionalInterface。

举例以Thread 为例说明很容易理解。

Runnable 接口是我们线程编程时常用的一个接口，就包含一个方法void run()，这个方法就是线程的运行逻辑。

按照以前的语法，我们新建线程一般要用到Runnable 的匿名类，如下：如果写多了，是不是很无聊，而基于Lambda 的写法则变得简洁明了，如下：注意Thread 的参数，Runnable 的匿名实现就通过一句就实现了出来，写成下面的更好理解当然Lambda 的目的不仅仅是写起来简洁，更高层次的目的等体会到了再总结。