CP计算物理常微分方程解实用.pptx
常微分方程描写的物理现象镭的衰变规律单摆的运动RLC振荡电路物理场计算第1页/共45页常微分方程第2页/共45页常微分方程数值解基本思想第3页/共45页8.1 Euler方法第4页/共45页8.1 Euler方法第5页/共45页8.1 Euler方法-梯形公式第6页/共45页8.2 改进Euler方法第7页/共45页8.2 改进Euler方法第8页/共45页12KC=0.01法R=10欧E=10伏第9页/共45页第10页/共45页function mainglobal E R C;E=10;R=10;C=0.01;Q0=E/C;h=0.01;t=0:h:1;QE(1)=Q0;QEG(1)=Q0;for i=2:length(t)QE(i)=QE(i-1)+h*f(QE(i-1);%欧拉法 k1=h*f(QEG(i-1);%改进欧拉法 k2=h*f(QEG(i-1)+k1);QEG(i)=QEG(i-1)+1/2*(k1+k2);endplot(t,QE,r);%欧拉法的曲线hold onplot(t,QEG,b);%改进欧拉法的曲线plot(t,Q0*exp(-1*t./(R*C),.);%理论曲线legend(欧拉法的曲线,改进欧拉法的曲线,理论曲线)function y=f(x)global E R C;y=-1*x./(R*C);第11页/共45页KC=10mFR=100欧E=10VL=10H 第12页/共45页第13页/共45页function main%欧拉法：二阶常微分方程global E R C L;E=10;R=100;C=0.01;L=10;Q0=0;I0=0;h=0.01;t=0:h:10;Q(1)=Q0;I(1)=I0;for i=2:length(t)Q(i)=Q(i-1)+h*I(i-1);%Q I(i)=I(i-1)+h*f(Q(i-1),I(i-1);%Iendplot(t,Q,r,t,I,m);%欧拉法的曲线function y=f(Q,I)global E R C L;y=(E-Q/C-I*R)/L;第14页/共45页（用改进的Euler法解）：第15页/共45页function main%改进欧拉法：二阶常微分方程global E R C L;E=10;R=100;C=0.01;L=10;Q0=0;I0=0;h=0.01;t=0:h:10;Q(1)=Q0;I(1)=I0;for i=2:length(t)%1、预报 Q(i)=Q(i-1)+h*I(i-1);I(i)=I(i-1)+h*f(Q(i-1),I(i-1);%2、计算Q k1=h*I(i-1);k2=h*I(i);Q(i)=Q(i-1)+1/2*(k1+k2);%3、计算I k1=h*f(Q(i-1),I(i-1);k2=h*f(Q(i),I(i);I(i)=I(i-1)+1/2*(k1+k2);endplot(t,Q,r,t,I,m);%改进欧拉法的曲线function y=f(Q,I)global E R C L;y=(E-Q/C-I*R)/L;第16页/共45页8.3 龙格-库塔（R-K）方法思想：取多点处斜率的加权平均为平均斜率，从而减小误差。四阶公式：公式推导见P37-44第17页/共45页例：用R-K方法解例题8.1.1第18页/共45页function main%龙格库塔：一阶常微分方程global E R C;E=10;R=10;C=0.01;h=0.01;t=0:h:1;Q(1)=E/C;for i=2:length(t)k1=h*f(Q(i-1);k2=h*f(Q(i-1)+k1/2);k3=h*f(Q(i-1)+k2/2);k4=h*f(Q(i-1)+k3);Q(i)=Q(i-1)+1/6*(k1+2*k2+2*k3+k4);endplot(t,Q,b);hold on%龙格库塔方法的曲线plot(t,Q(1)*exp(-1*t./(R*C),.);%理论曲线function y=f(x)global E R C;y=-1*x./(R*C);第19页/共45页例：用R-K方法解例题8.1.2和第20页/共45页function main%龙格库塔：二阶常微分方程global E R C L;E=10;R=100;C=0.01;L=10;Q0=0;I0=0;h=0.01;t=0:h:10;Q(1)=Q0;I(1)=I0;for i=2:length(t)k1=h*I(i-1);m1=h*f(Q(i-1),I(i-1);k2=h*(I(i-1)+m1/2);m2=h*f(Q(i-1)+k1/2,I(i-1)+m1/2);k3=h*(I(i-1)+m2/2);m3=h*f(Q(i-1)+k2/2,I(i-1)+m2/2);k4=h*(I(i-1)+m3);m4=h*f(Q(i-1)+k3,I(i-1)+m3);Q(i)=Q(i-1)+1/6*(k1+2*k2+2*k3+k4);I(i)=I(i-1)+1/6*(m1+2*m2+2*m3+m4);endplot(t,Q,r,t,I,m);%龙格库塔的曲线function y=f(Q,I)global E R C L;y=(E-Q/C-I*R)/L;第21页/共45页练习：分别用Euler法、改进Euler法和四阶R-K法求解阻尼振动方程：已知质量m=10,倔强系数k=10,阻尼系数c=2,初始速度v=0，初始位置x=10.第22页/共45页误差概述第23页/共45页误差概述第24页/共45页误差概述第25页/共45页误差概述第26页/共45页8.1.3 数值稳定性分析第27页/共45页数值稳定性分析定义8.1.3 若某数值算法的绝对稳定性区域包含h平面上的左半平面Re(h)0，则称该方法是A稳定的。隐式Euler法是A稳定的。第28页/共45页8.2 Runge-Kutta方法第29页/共45页Runge-Kutta方法第30页/共45页Runge-Kutta方法第31页/共45页第32页/共45页Runge-Kutta方法第33页/共45页第34页/共45页6.2.2 四阶Runge-Kutta方法第35页/共45页四阶Runge-Kutta方法第36页/共45页6.2.3 R-K法的稳定性第37页/共45页R-K法的稳定性第38页/共45页R-K法的稳定性第39页/共45页6.2.5 隐式R-K法第40页/共45页隐式R-K法第41页/共45页隐式R-K法第42页/共45页隐式R-K法第43页/共45页隐式R-K法第44页/共45页感谢您的欣赏！第45页/共45页
