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

一阶常微分方程的数值求解; 考虑如下的一维经典初值问题; 具体步骤：;例1：用 Euler 法求解如下初值问题;clc;clear; h=0.1; a=0;b=2; x=a:h:b; y(1)=1; for i=1:length(x)-1 y(i+1)=y(i)+h*0.5*y(i); end plot(x,y,r+,x,exp(0.5*x),k-); xlabel(Variable x); ylabel(Variable y); ;数值解与真解如下图;例2：用 Euler 法求解如下初值问题;clc;clear; h=0.1; a=0;b=2; x=a:h:b; y(1)=1; for i=1:length(x)-1 y(i+1)=y(i)+h*(y(i)+2*x(i)/(y(i)^2)); end plot(x,y,r+,x,(5/3*(exp(3*x))-2*x-2/3).^(1/3),k-); xlabel(Variable x); ylabel(Variable y); ;数值解与真解如下图; 为了减小误差，可采用以下方法：; 常用的是经典的 四阶R-K方法;例3：利用四阶R-K方法求解例1与例2，并与Euler方法的数值解进行比较。;clc;clear; h=0.1; a=0;b=2; x=a:h:b; Euler_y(1)=1; % Euler方法的初值 RK_y(1)=1; % R-K方法的初值 for i=1:length(x)-1 Euler_y(i+1)=Euler_y(i)+h*0.5*Euler_y(i); %Euler方法 L1=0.5*RK_y(i); % L1 L2=0.5*(RK_y(i)+0.5*h*L1); % L2 L3=0.5*(RK_y(i)+0.5*h*L2); % L3 L4=0.5*(RK_y(i)+h*L3); % L4 RK_y(i+1)=RK_y(i)+1/6*h*(L1+2*L2+2*L3+L4); %R-K方法 end plot(x,Euler_y,r+,x,exp(0.5*x),k-,x,RK_y,b*); xlabel(Variable x); ylabel(Variable y); ;当 h=0.1，即 n=20 时，求解例2的Matlab 源程序见 RK_example2.m, 数值结果如下图;function z=rightf(x,y) z=y+2*x/(y^2);;Matlab函数数值求解;Matlab提供的ODE求解器;作业