药物计算分析导论章.pptx

会计学1药物计算分析导论章药物计算分析导论章2003年This LectureThis Lecturen nMore accurate schemesMore accurate schemesn nMore complicated ODEsMore complicated ODEsn nVariable time step and embedded methods used to make sure Variable time step and embedded methods used to make sure errors are within a tolerance.errors are within a tolerance.第1页/共46页2003年Adams-Bashforth SchemesAdams-Bashforth Schemesn nIn the forward Euler scheme we only used the value of the right hand side at the In the forward Euler scheme we only used the value of the right hand side at the previous time step.previous time step.n ni.e.we only used a linear approximation to the time derivativei.e.we only used a linear approximation to the time derivative第2页/共46页2003年AB SchemesAB Schemesn nIdea:we set Idea:we set n nwhere If interpolates fwhere If interpolates fn n,f,fn-1n-1,f,fn-2n-2,.,f,.,fn+1-Nstagesn+1-Nstages n ni.e.:i.e.:第3页/共46页2003年f0f1f2f3We interpolate the function through thefirst 4 points.Then we integrate under the curve between t=3 and t=4.第4页/共46页2003年AB SchemesAB Schemes Essentially we use Essentially we use interpolation and a interpolation and a Newton-Cotes quadrature Newton-Cotes quadrature formula to formulate:formula to formulate:第5页/共46页2003年Runge-Kutta SchemesRunge-Kutta Schemesn nSee van Loan for derivation of Runge-Kutta2 and Runge-See van Loan for derivation of Runge-Kutta2 and Runge-Kutta4.Kutta4.n nThe following(simple)scheme due to Jameson,Schmidt The following(simple)scheme due to Jameson,Schmidt and Turkel(1981):and Turkel(1981):第6页/共46页2003年Runge-Kutta SchemesRunge-Kutta Schemesn nBeware,it only works when Beware,it only works when f f is a function of is a function of y y and and notnot t tn nhere here s s is the order of the scheme.is the order of the scheme.第7页/共46页2003年Error EstimateError Estimaten nMatlab has a number of time integrators built in.Matlab has a number of time integrators built in.n node23ode23n node45ode45n nand others.and others.第8页/共46页2003年ode23ode23n nFor n=1:#timestepsFor n=1:#timestepsn node23 uses two estimates for yode23 uses two estimates for yn+1n+1.n nA 2A 2ndnd order RK scheme and a 3 order RK scheme and a 3rdrd order RK scheme are used to build two order RK scheme are used to build two guesses for yguesses for yn+1n+1.n nIf the difference between these two estimates are within a tolerance ode23 If the difference between these two estimates are within a tolerance ode23 progresses on to calculating yprogresses on to calculating yn+2n+2n nIf the difference is greater than the specified tolerance,ode23 reduces the dt If the difference is greater than the specified tolerance,ode23 reduces the dt and tries again.It repeats until the difference is lower than the tolerance.and tries again.It repeats until the difference is lower than the tolerance.n nEndEnd第9页/共46页2003年T,Y=ode23(odefun,tspan,yT,Y=ode23(odefun,tspan,y0 0)nwith tspan=t0 tf integrates the system of differential equations from time t0 to tf with initial conditions y0.nY0 A vector of initial conditions.nOdefun A function that evaluates the right-hand side of the differential equations.第10页/共46页2003年Planets Example Using ode23Planets Example Using ode23n nIdea:replace our home grown Euler Forward scheme with MatlabIdea:replace our home grown Euler Forward scheme with Matlab s ode23 s ode23 scheme in the scheme in the planets1.mplanets1.m script.script.第11页/共46页2003年Parameters forthe planets asbefore第12页/共46页2003年Initial velocitiesGatherX,Y,VX,VYinto one vector第13页/共46页2003年 T,CoordVel=ode23(forcing,TimeLimits,CoordVel(:);Call to Matlab ode23 functionfunction whichcalculates f(y)0 FinalTime Vector holdsX,Y,VX,VYWorks in Matlab v6.1.0 at least第14页/共46页2003年Call to ode23Find out the time at each time stepExtract final coordinates and velocitiesPlot planet orbitsd=size(X)returns the sizes of each dimension of array X in a vector d with ndims(X)elements.第15页/共46页2003年The routine whichcalculates the forcing function.第16页/共46页2003年Team ExerciseTeam Exercisen nGrab Grab planets2.mplanets2.m and and forcing.mforcing.mn nRun the scriptRun the scriptn nUse the Tsteps vector to find out the time step for each Use the Tsteps vector to find out the time step for each integration stage.Plot a graph showing the time step(dt)integration stage.Plot a graph showing the time step(dt)at each time step.at each time step.n nUse Use helphelp to find out how to change the tolerance used by to find out how to change the tolerance used by ode23ode23 (hint you will need to use (hint you will need to use odesetodeset)n nRerun the simulation with a tolerance of 0.1 Rerun the simulation with a tolerance of 0.1 第17页/共46页2003年第18页/共46页2003年Application:One-Dimensional Electrostatic Application:One-Dimensional Electrostatic MotionMotion第19页/共46页2003年Charge RepulsionCharge Repulsionn nNow we will consider the case of charged particles with the same sign Now we will consider the case of charged particles with the same sign chargechargen nInstead of attracting each other,the charges repel each other.Instead of attracting each other,the charges repel each other.第20页/共46页2003年Particle AccelerationsParticle Accelerationsa1a2Each particle has a vector accelerationdirectly away from the other particle第21页/共46页2003年Team ProjectTeam ProjectQ1)Modify the Q1)Modify the planets2.mplanets2.m and and forcing.mforcing.m to simulate the following:to simulate the following:n nThere are N electrically charged particles There are N electrically charged particles confined to move in the x-confined to move in the x-direction onlydirection onlyn nDistribute the charges initially at equispaced points in-1,1Distribute the charges initially at equispaced points in-1,1n nThe equations of motion of the charges are:The equations of motion of the charges are:第22页/共46页2003年Team ProjectTeam ProjectQ1cont)Include the following in one project write up per team:Q1cont)Include the following in one project write up per team:n nA scatter plot,with x as the horizontal axis and t as the vertical axis,showing A scatter plot,with x as the horizontal axis and t as the vertical axis,showing the paths of all the charged particles using ode23the paths of all the charged particles using ode23n nReplace ode23 with ode45 and rerunReplace ode23 with ode45 and rerunn nPlot the ode23 and ode45(t,x)paths on the same graph.Plot the ode23 and ode45(t,x)paths on the same graph.n nPlot the ode23 and ode45 time step sizes on the same graphPlot the ode23 and ode45 time step sizes on the same graphn nNames of team membersNames of team members第23页/共46页2003年Chap.12 Chap.12 Electrostatics-Colliding Disks ProjectElectrostatics-Colliding Disks Project第24页/共46页2003年Big Team ProjectBig Team Projectn nOk Ok we have been warming up with some small projects we have been warming up with some small projectsn nNow for something a little more involvedNow for something a little more involved第25页/共46页2003年The PhysicsThe Physicsn nThe domain of interest is a two-dimensional periodic The domain of interest is a two-dimensional periodic box of side length Lbox of side length Ln nIt is populated with N circles of mass M each and It is populated with N circles of mass M each and radius Rradius R第26页/共46页2003年xyx=0y=0 x=Ly=LDomain第27页/共46页2003年xyx=0y=0 x=Ly=LParticlesR第28页/共46页2003年Defining A Set of Motion RulesDefining A Set of Motion RulesRule 1Rule 1The disks are initially randomly distributedThe disks have random initial velocity(with the absolute value of each component of velocity bounded by 1)Rule 2Rule 2Effectively Newtons first:all particles do not accelerate until they collide第29页/共46页2003年Particles CollideParticles Collide21Before21After第30页/共46页2003年Rule 4Rule 4n nMomentum(Momentum(动量动量)is conserved in a collision)is conserved in a collisionRule 5Rule 5Angular momentum(动量矩)is conserved in the collisionRule 6Rule 6Total kinetic energy(动能)is conserved The total kinetic energy of the particles is the same before and after the collision.第31页/共46页2003年Model SimplificationModel Simplificationn nWe are going to use Rules 4,5,and 6 to determine the change in We are going to use Rules 4,5,and 6 to determine the change in velocity of two disks when they collide.velocity of two disks when they collide.n nWe can simplify this problem by considering a coordinate frame We can simplify this problem by considering a coordinate frame which is:which is:n n co-moving with the center of one of the disks(say disk1)at co-moving with the center of one of the disks(say disk1)at its originits originn n with x-axis aligned in direction of the segment connecting the with x-axis aligned in direction of the segment connecting the two disk centerstwo disk centers第32页/共46页2003年In Co-Moving FrameIn Co-Moving Frame21Before第33页/共46页2003年Math Version of 4Math Version of 4(in co-moving frame)(in co-moving frame)n nConservation of Conservation of linear momentumlinear momentum:n nWhich reduces to:Which reduces to:第34页/共46页2003年Math Version of 5Math Version of 5(in co-moving frame)(in co-moving frame)n nConservation of Conservation of angular angular momentum:momentum:n nAnd since the center of two can not be at the origin And since the center of two can not be at the origin this implies:this implies:第35页/共46页2003年Math Version of 6Math Version of 6(in co-moving frame)(in co-moving frame)n nConservation of Conservation of kinetic energykinetic energy:n nReduces to:Reduces to:第36页/共46页2003年Summary of 4,5,6Summary of 4,5,6456第37页/共46页2003年SolutionsSolutionsCollisionNo collision第38页/共46页2003年In Co-Moving FrameIn Co-Moving Frame21Before21After第39页/共46页2003年Reverting To Non-moving FrameReverting To Non-moving Frame The rules essentially say that if two equal The rules essentially say that if two equalmass disks collide then:mass disks collide then:1)1)they exchange velocity in their center-to-center vector they exchange velocity in their center-to-center vector 2)2)they retain their own velocity in the direction of the tangent to they retain their own velocity in the direction of the tangent to their center-to-center vectortheir center-to-center vector第40页/共46页2003年General FormulaGeneral Formula*第41页/共46页2003年Summary of Collision FormulaSummary of Collision Formulan nGiven the relative position of two disks before the collision and their velocities we can determine the post-collision velocities第42页/共46页2003年Change invelocity ofdisk 1Exchange of linear andangular momentum第43页/共46页2003年Project DescriptionProject Descriptionn nUsing Euler forward as time steppingUsing Euler forward as time steppingn nInitiate the random initial positions and velocities of the disksInitiate the random initial positions and velocities of the disksn nMove the particles with time step dtMove the particles with time step dtn nEach time step check to see if any disk intersects any other diskEach time step check to see if any disk intersects any other diskn nif two disks intersect use the collision formula to determine their new velocitiesif two disks intersect use the collision formula to determine their new velocitiesn nIncrement(Increment(增大增大)the disk positions with dt multiplying their velocitiesthe disk positions with dt multiplying their velocities第44页/共46页2003年Project Description contProject Description cont Use mod on the updated positions to make sure that the disks stay within the box Repeat until quite a few collisions show up Plot the locations of the particles in a two-dimensional plot use hold so that the time history of the disks shows upmod(x,y)is X X mod Y Y 第45页/共46页