数学实验2015_05.pdfVIP

数学实验 张云飞 计算数学教研室 Tel: 微分方程数值解 算法； 欧拉与龙格-库塔方法 模型与实验 算法 当常微分方程能解析求解时，可利用Matlab符号 工具箱中的功能找到精确解。见下例求解方程。 syms x y % 定义符号变量 diff_equ = ‘D2y +2*Dy -y =0’; % D2y表示y = dsolve(diff_equ,’x’) % 定义x为自变量 % y = % c1*exp((2^(1/2)-1)*x)+c2*exp(-(2^(1/2)+1)*x) % 初始条件y(0) = 0, y’(0) = 1， 按如下方式调用 y = dsolve(diff_equ,’y(0)=0’,’Dy(0)=1’,’x’) % y = % 1/4*2^(1/2)*exp((2^(1/2)-1)*x)-1/4*2^(1/2)*exp(- (1+2^(1/2))*x) % 画出函数y = y(x) 的图形 ezplot(y,[-2,2]) 方程无解析解时，可获得数值解。 下例说明Matlab求数值解方法及注意问题。 [例1] 求解范德堡（vander pol ）方程 d 2 x 2 dx −µ(1− ) +x 0 x dt 2 dt 求解高阶方程，必须等价地变换为一阶微分方程组， 本例，通过定义两个新的变量，实现这一变换 y 1 x , y 2 dx dt 则令 dy1 dt y 2 2 dy2 dt µ(1−y 1 ) ∗y 2 −y 1 求解程序分两部分，第一部分为待求解的方程。存 盘的文件名为待求解方程的函数名.m，第二部分 为求解主程序。 function yprime = vdpol(t,y ) yprime(1) = y(2); mu = 2; yprime(2) = mu*(1-y(1)^2)*y(2)-y(1); yprime = [yprime(1);yprime(2)]; 说明 函数yprime = vdpol(t,y ) ，t为自变量，y 的形式 取决于求解方程的阶数，y = [y (1), y (2)]，y (1)为解向 量，y (2)为导数向量。yprime(1) = y ’(1)，yprime(2) = y ”(1)，函数返回vander pol方程的导数列向量。因为 所求结果为方程数值解，所以各向量维数只有在主 程序求解时定下精度后才能确定。 主程序main1.m 。 clear functions % 清除前边已编译的存在于内存中的废弃程序 [t,y ] = ode23(‘vdpol’,[0,30],[1,0]); y1 = y(:,1); % 解曲线。 y2 = y(:,2); % 解曲线的导数。 plot(t,y 1,t,y 2,’--’) 说明 龙格_库塔的2阶与4阶改进型求解公式： [t,x] = ode23(‘f ’, ts, x0, options) [t,x] = ode45(‘f ’, ts, x0, options) 其中ts可由系统精度要求自动设定，亦可由使用者 自己确定。 (1) 当ts = [t0, t1, …, tf]，指定时刻t0, t1, …, tf给出 输出，当ts = t0: k: tf ，区间[t0, tf]等分点上输出。 (2) 若ts = [t0, tf]，options决定自变量t 的维数，t 中 的时间点非等间隔，保证所需的相对精度，积分算 法改变了步长。误差参数options缺省时设定相对误 差10-3 ，绝对误差10-6 ，若自行设定，用如下语句： options = o