CP计算物理常微分方程解实用实用教案

1、常微分方程(wi fn fn chn)描写的物理现象 镭的衰变规律 单摆的运动 RLC振荡电路 物理(wl)场计算第1页/共45页第一页，共46页。常微分方程(wi fn fn chn)第2页/共45页第二页，共46页。常微分方程(wi fn fn chn)数值解基本思想第3页/共45页第三页，共46页。8.1 Euler方法(fngf)第4页/共45页第四页，共46页。8.1 Euler方法(fngf)第5页/共45页第五页，共46页。8.1 Euler方法(fngf)-梯形公式第6页/共45页第六页，共46页。8.2 改进(gijn)Euler方法第7页/共45页第七页，共46页。8.2

2、改进(gijn)Euler方法第8页/共45页第八页，共46页。12KC=0.01法R=10欧E=10伏第9页/共45页第九页，共46页。第10页/共45页第十页，共46页。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(

3、t,QE,r);%欧拉法的曲线(qxin)hold onplot(t,QEG,b);%改进欧拉法的曲线(qxin)plot(t,Q0*exp(-1*t./(R*C),.);%理论曲线(qxin)legend(欧拉法的曲线(qxin),改进欧拉法的曲线(qxin), 理论曲线(qxin)function y=f(x)global E R C;y=-1*x./(R*C);第11页/共45页第十一页，共46页。KC=10mFR=100欧E=10VL=10H 第12页/共45页第十二页，共46页。第13页/共45页第十三页，共46页。function main%欧拉法：二阶常微分方程(wi fn fn

4、 chn)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页第十四页，共46页。（用改进(gijn)的Euler法解）：第15页/共45页第十五页，共46页。functio

5、n main%改进欧拉法：二阶常微分方程(wi fn fn chn)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

6、(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页第十六页，共46页。8.3 龙格-库塔（R-K）方法(fngf) 思想(sxing)：取多点处斜率的加权平均为平均斜率，从而减小误差。 四阶公式： 公式(gngsh)推导见P37-44第17页/共45页第十七页，共46页。 例：用R-K方法(fngf)解例题8.1.1第18页/共45页第十八页，共46页。function main%龙格库塔：一阶常微分方程global E R C;E=1

7、0;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%龙格库塔方法的曲线(qxin)plot(t,Q(1)*exp(-1*t./(R*C),.);%理论曲线(qxin)function y=f(x)global E R C;y=-1*x./(R*C);第19页/共45页第十九

8、页，共46页。 例：用R-K方法(fngf)解例题8.1.2和第20页/共45页第二十页，共46页。function main%龙格库塔：二阶常微分方程(wi fn fn chn)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

9、*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页第二十一页，共46页。练习(linx)： 分别用Euler法、改进Euler法和四阶R-K法求解阻尼振动方程： 已知质量m=10,倔强系数

10、k=10,阻尼系数c=2,初始(ch sh)速度v=0，初始(ch sh)位置x=10.第22页/共45页第二十二页，共46页。误差(wch)概述第23页/共45页第二十三页，共46页。误差(wch)概述第24页/共45页第二十四页，共46页。误差(wch)概述第25页/共45页第二十五页，共46页。误差(wch)概述第26页/共45页第二十六页，共46页。8.1.3 数值(shz)稳定性分析第27页/共45页第二十七页，共46页。数值(shz)稳定性分析 定义8.1.3 若某数值算法的绝对稳定性区域(qy)包含h平面上的左半平面Re(h)0，则称该方法是A稳定的。 隐式Euler法是A稳定的

11、。第28页/共45页第二十八页，共46页。8.2 Runge-Kutta方法(fngf)第29页/共45页第二十九页，共46页。Runge-Kutta方法(fngf)第30页/共45页第三十页，共46页。Runge-Kutta方法(fngf)第31页/共45页第三十一页，共46页。第32页/共45页第三十二页，共46页。Runge-Kutta方法(fngf)第33页/共45页第三十三页，共46页。第34页/共45页第三十四页，共46页。6.2.2 四阶Runge-Kutta方法(fngf)第35页/共45页第三十五页，共46页。四阶Runge-Kutta方法(fngf)第36页/共45页第三十

12、六页，共46页。6.2.3 R-K法的稳定性第37页/共45页第三十七页，共46页。R-K法的稳定性第38页/共45页第三十八页，共46页。R-K法的稳定性第39页/共45页第三十九页，共46页。6.2.5 隐式R-K法第40页/共45页第四十页，共46页。隐式R-K法第41页/共45页第四十一页，共46页。隐式R-K法第42页/共45页第四十二页，共46页。隐式R-K法第43页/共45页第四十三页，共46页。隐式R-K法第44页/共45页第四十四页，共46页。感谢您的欣赏(xnshng)！第45页/共45页第四十五页，共46页。NoImage内容(nirng)总结常微分方程描写的物理现象。第1页