常微分方程实验报告_闵安游_U201310085

⑵、3 步 4 阶 Adams 方法： ①、画图程序：
w = exp(i*linspace(0,2*pi)); z = 24*(w.^3 - w.^2)./(1 - 5*w + 19*w.^2 + 9*w.^3); plot(z)
②、绝对稳定性判断程序：
9 3 19 2 5 h (1 h ) (1 h) h 0 此时的特征方程是： 24 24 24 24
fprintf('误差为: %10.10f\n',m)
测试程序：
f = @(t,y)-2*y + 2*t + 1; F = @(t)exp(-2*t) + t; Threesteps_Adams1(0,2,f,F,0.1)
第 2 题： ⑴、2 步 3 阶 Adams 方法： ①、画图程序：
w = exp(i*linspace(0,2*pi)); z = 12*w.*(w - 1)./(5*w.*w + 8*w -1); plot(z)
t = a + i*h; y2 = y1 + h/2*(3*f1 - f0); plot(t,y2,'.'); f2 = f(t,y2); y1 = y2; f0 = f1; f1 = f2; end y = linspace(a,b); plot(y,F(y),'r'); hold off m = abs(y1-F(b)); fprintf('数值解: %10.10f\n',y1) fprintf('误差为:%10.10f\n',m)
' u2 24u1 51u2 9 cos t 1/ 3sin t
u1 (0) 4 / 3, u2 2 / 3
列出当步长 h 0.2, 0.1, 0.05 时计算到终点 t 2 的数值解及其误差，其 中初始值直接采用准确值，该问题的真解是
u1 (t ) 2e 3t e 39t 1/ 3cos t u2 (t ) e 3t 2e 39t 1/ 3cos t
测试程序：
f = @(t,y)-2*y + 2*t + 1; F = @(t)exp(-2*t) + t; Twosteps_Adams1(0,2,f,F,0.2)
⑵、3 步 3 阶显式 Adams 方法： ①、初值由真解给定： 主程序：
function Threesteps_Adams(a,b,f,F,h) %初值由真解给出 N = (b - a)/h; y0 = F(a); y1 = F(a + h); y2 = F(a + 2*h); plot(a,y0,'.') hold on
要求对所获得结果进行分析，并附两个算法的 MATLAB 程序。 4、应用显示中点方法和 4 级 4 阶 Gill 方法求解初值问题
y ' y t , t [0,1], y (0) 3
列出当步长 h 0.2 和 0.1 时计算到终点 t 1 的数值解及其误差，并附
t 算法的 MATLAB 程序。 （问题的真解为 y (t ) 2e t 1 ）
第 3 题： ⑴、2 步 3 阶隐式 Adams 方法： ①、主程序：
function x = TwoStep3jie(h) %定义给定具体的 h 的函数 a = 0;b = 2; n = (b - a)/h; I = eye(2); A = [9 24;-24 -51];% 关于向量 u 的向量矩阵 B = [5 -1/3;-9 1/3];% 关于时间 t 的向量矩阵 u0 = F2(0); u1 = F2(a + h); subplot(1,2,1) plot(a,u0(1,1),'.') hold on plot(a + h,u1(1,1),'.') hold on subplot(1,2,2) plot(a,u0(2,1),'.') hold on
测试程序：
f = @(t,y)-2*y + 2*t + 1; F = @(t)exp(-2*t) + t; Threesteps_Adam程序：
function Threesteps_Adams1(a,b,f,F,h) %初值由显示欧拉给出 N = (b - a)/h; y0 = F(a); f0 = f(a,y0); y1 = y0 + h*f0; f1 = f(a + h,y1); y2 = y1 + h*f1; f2 = f(a + 2*h,y2); plot(a,y0,'.') hold on plot(a+h,y1,'.') hold on plot(a + 2*h,y2,'.') hold on for i = 3:N t = a + i*h; y3 = y2 + h/12*(23*f2 - 16*f1 + 5*f0); plot(t,y3,'.') f3=f(t,y3); y2=y3; f0=f1; f1=f2; f2=f3; end y = linspace(a,b); plot(y,F(y),'r'); hold off m = abs(y2-F(b)); fprintf('数值解: %10.10f\n',y2)

所以程序为：
function test34(h) p = [1-9/24*h 1+19/24*h h*5/24 -h/24]; m = roots(p); fprintf('多项式根为: %5.5f\n',m(1,1)) fprintf(' fprintf(' %5.5f\n',m(2,1)) %5.5f\n',m(3,1))
测试程序：
f = @(t,y)-2*y + 2*t + 1; F = @(t)exp(-2*t) + t; Twosteps_Adams(0,2,f,F,0.2)
②、初值由显示欧拉给定： 主程序：
function Twosteps_Adams1(a,b,f,F,h) N = (b - a)/h; y0 = F(a); f0 = f(a,y0); y1 = y0 + h*f0; f1 = f(a + h,y1); plot(a,y0,'.') hold on plot(a+h,y1,'.') hold on for i = 2 :N
②、绝对稳定性判断程序：
5 2 h 此时的特征方程是： (1 h) 2 (1 h) 0 12 3 12

所以程序为：
function test23(h) p = [1-5/12*h 1+2/3*h h/12]; m = roots(p); fprintf('多项式根为: %5.5f\n',m(1,1)) fprintf(' %5.5f\n',m(2,1))
plot(a+h,y1,'.') hold on plot(a + 2*h,y2,'.') hold on f0 = f(a,y0); f1 = f(a + h,y1); f2 = f(a + 2*h,y2); for i = 3:N t = a + i*h; y3 = y2 + h/12*(23*f2 - 16*f1 + 5*f0); plot(t,y3,'.') f3=f(t,y3); y2=y3; f0=f1; f1=f2; f2=f3; end y = linspace(a,b); plot(y,F(y),'r'); hold off m = abs(y2-F(b)); fprintf('数值解: %10.10f\n',y2) fprintf('误差为: %10.10f\n',m)
2. 实验题目
1、分别应用 2 步 2 阶与 3 步 3 阶显式 Adams 方法求解初值问题
y ' 2 y 2t 1, t [0, 2], y (0) 1
列出当步长 h 0.2 和 0.1 时计算到终点 t 2 的数值解及其误差，其中 的初值考虑两种方案，一种是使用显式 Euler 方法，另一种是直接 采用准确值(其真解为 y e
华中科技大学
本科实验报告
课程名称：常微分方程数值解 姓 学 名：闵安游 号：U201310085
专业班级: 信计 1301 班 指导老师：黄乘明 实验题目：常微分方程数值解上机 日 期：2015.12.25 实验成绩：
1. 实验目的
通过上机实验，理解利用计算机迭代求解线性方程组的整个过 程，加深对所学计算方法的理论及算法特点的理解。
5、应用 3 级 3 阶 Heun 方法求解初值问题
y ' 2 y 2t 1, t [0, 2], y (0) 1
列出当步长 h 0.2 和 0.1 时计算到终点 t 2 的数值解及其误差，并附 算法的 MATLAB 程序.(该问题的真解为 y e 2t t ，其中方法所需要 的初值直接使用准确值)
plot(a + h,u1(2,1),'.') hold on for i = 0:n - 2; t = a + i*h; u2 = (I-(5*h/12)*A)^(-1)*(u1 + (5*h/12)*F3(t+2*h) + (h/12)*(8*F1(u1,t+h) -F1(u0,t)));% 迭代项 subplot(1,2,1) plot(t + 2*h,u2(1,1),'.') hold on subplot(1,2,2) plot(t + 2*h,u2(2,1),'.') u0 = u1; u1 = u2; end u2; y1=abs(F2(2) - u2); fprintf('u1 数值解: %10.10f\n',u2(1,1)) fprintf('u2 数值解: %10.10f\n',u2(2,1)) fprintf('u1 误差为:%10.10f\n',y1(1,1)) fprintf('u2 误差为:%10.10f\n',y1(2,1)) y = linspace(a,b); subplot(1,2,1) plot(y,2*exp(-3*y)-exp(-39*y)+1/3*cos(y),'r')% u1 准确值的图像 （红色的线表示） hold on subplot(1,2,2) plot(y,-exp(-3*y)+2*exp(-39*y)-1/3*cos(y),'r')% u2 准确值的图像 （红色的线表示）
2 t
t )。要求对所获得结果进行分析，并
附两个算法的 MATLAB 程序。 2、画出 2 步 3 阶和 3 步 4 阶隐式 Adams 方法的边界轨迹，并确定 方法的绝对稳定区域。 3、应用 2 步 3 阶和 3 步 4 阶隐式 Adams 方法求解初值问题
u1' 9u1 24u2 5cos t 1/ 3sin t
3. 程序文本
第 1 题: ⑴、2 步 2 阶 Adams 方法： ①、初值由真解给定： 主程序：
function Twosteps_Adams(a,b,f,F,h) % a,b 为区间端点，F 为真解函 数，h 为步长。 N = (b - a)/h; y0 = F(a); y1 = F(a + h); plot(a,y0,'.') hold on plot(a+h,y1,'.') hold on f0 = f(a,y0); f1 = f(a + h,y1); for i = 2 :N t = a + i*h; y2 = y1 + h/2*(3*f1 - f0);
②、调用程序：
function z = F1(u,t) % 表示 f(t,u)这个向量函数 A = [9 24;-24 -51]; B = [5 -1/3;-9 1/3]; z = A*u + B*[cos(t);sin(t)];
plot(t,y2,'.'); f2 = f(t,y2); y1 = y2; f0 = f1; f1 = f2; end y = linspace(a,b); plot(y,F(y),'r');%解析函数图像 hold off m = abs(y2-F(b));% 误差 fprintf('数值解: %10.10f\n',y1) fprintf('误差为: %10.10f\n',m)