一阶常微分方程的数值求解

一阶常微分方程的数值求解,一.教学要求,掌握Euler方法数值求解一阶常微分初值问题，并能利用MATLAB软件进行数值计算。,二.教学过程,考虑如下的一维经典初值问题,基本思想：用差商代替微商,根据Talyor公式，y(x)在点xk处有,具体步骤：,等距剖分：,步长：,分割求解区间,差商代替微商,得方程组：,分割求解区间，差商代替微商，解代数方程,为分割点,k=0,1,2,.,n-1,yk是y(xk)的近似,例1：用Euler法求解如下初值问题,取步长h=(2-0)/n=2/n，得差分方程,当h=0.1，即n=20时，Matlab源程序见Euler_example1.m,解：,clc;clear;h=0.1;a=0;b=2;x=a:h:b;y(1)=1;fori=1:length(x)-1y(i+1)=y(i)+h*0.5*y(i);endplot(x,y,r+,x,exp(0.5*x),k-);xlabel(Variablex);ylabel(Variabley);,Euler_example1.m,数值解与真解如下图,例2：用Euler法求解如下初值问题,取步长h=(2-0)/n=2/n，得差分方程,当h=0.1，即n=20时，Matlab源程序见Euler_example2.m,解：,clc;clear;h=0.1;a=0;b=2;x=a:h:b;y(1)=1;fori=1:length(x)-1y(i+1)=y(i)+h*(y(i)+2*x(i)/(y(i)2);endplot(x,y,r+,x,(5/3*(exp(3*x)-2*x-2/3).(1/3),k-);xlabel(Variablex);ylabel(Variabley);,Euler_example2.m,数值解与真解如下图,为了减小误差，可采用以下方法：,让步长h取得更小一些；如取h=0.01,h=0.001,改用具有较高精度的数值方法：,龙格-库塔方法,如Runge-Kutta(龙格-库塔)方法,是一类求解常微分方程的数值方法,有多种不同的迭代格式,龙格-库塔方法的基本思想,常用的是经典的四阶R-K方法,其中,例3：利用四阶R-K方法求解例1与例2，并与Euler方法的数值解进行比较。,当h=0.1，即n=20时，求解例1的Matlab源程序见RK_example1.m,数值结果如下图,clc;clear;h=0.1;a=0;b=2;x=a:h:b;Euler_y(1)=1;%Euler方法的初值RK_y(1)=1;%R-K方法的初值fori=1:length(x)-1Euler_y(i+1)=Euler_y(i)+h*0.5*Euler_y(i);%Euler方法L1=0.5*RK_y(i);%L1L2=0.5*(RK_y(i)+0.5*h*L1);%L2L3=0.5*(RK_y(i)+0.5*h*L2);%L3L4=0.5*(RK_y(i)+h*L3);%L4RK_y(i+1)=RK_y(i)+1/6*h*(L1+2*L2+2*L3+L4);%R-K方法endplot(x,Euler_y,r+,x,exp(0.5*x),k-,x,RK_y,b*);xlabel(Variablex);ylabel(Variabley);,RK_example1.m,当h=0.1，即n=20时，求解例2的Matlab源程序见RK_example2.m,数值结果如下图,functionz=rightf(x,y)z=y+2*x/(y2);,clc;clear;h=0.1;a=0;b=2;x=a:h:b;Euler_y(1)=1;%Euler方法的初值RK_y(1)=1;%R-K方法的初值fori=1:length(x)-1Euler_y(i+1)=Euler_y(i)+h*(Euler_y(i)+2*x(i)/(Euler_y(i)2);%Euler方法L1=rightf(x(i),RK_y(i);%L1L2=rightf(x(i)+0.5*h,RK_y(i)+0.5*h*L1);%L2L3=rightf(x(i)+0.5*h,RK_y(i)+0.5*h*L2);%L3L4=rightf(x(i)+h,RK_y(i)+h*L3);%L4RK_y(i+1)=RK_y(i)+1/6*h*(L1+2*L2+2*L3+L4);%R-K方法endplot(x,Euler_y,r+,x,(5/3*(exp(3*x)-2*x-2/3).(1/3),k-,x,RK_y,b*);xlabel(Variablex);ylabel(Variabley);,rightf.m,RK_example2.m,Matlab函数数值求解,T,Y=solver(odefun,tspan,y0),其中y0为初值条件，tspan为求解区间；Matlab在数值求解时自动对求解区间进行分割，T(向量)中返回的是分割点的值(自变量)，Y(向量)中返回的是解函数在这些分割点上的函数值。solver为Matlab的ODE求解器（可以是ode45、ode23、ode113、ode15s、ode23s、ode23t、ode23tb）,没有一种算法可以有效地解决所有的O