104753150620 姚鸿泰 关于SI模型的线性多步法求解讨论.docx
- 文档编号:27624963
- 上传时间:2023-07-03
- 格式:DOCX
- 页数:11
- 大小:132.14KB
104753150620 姚鸿泰 关于SI模型的线性多步法求解讨论.docx
《104753150620 姚鸿泰 关于SI模型的线性多步法求解讨论.docx》由会员分享,可在线阅读,更多相关《104753150620 姚鸿泰 关于SI模型的线性多步法求解讨论.docx(11页珍藏版)》请在冰豆网上搜索。
104753150620姚鸿泰关于SI模型的线性多步法求解讨论
高等数值分析课程论文
关于SI模型的线性多步法求解讨论
姓名:
姚鸿泰
学号:
104753150620
学院:
数学与统计学院
专业:
概率论与数理统计
时间:
2016年6月
关于SI模型的线性多步法求解讨论
姚鸿泰概率论与数理统计104753150620
摘要:
本文首先介绍了SI模型——易感染者与已感染者模型,并对该模型建模,得到非线性常微分方程.之后运用线性多步法求解该非线性常微分方程,并在最终给出误差图.
关键词:
SI模型;线性多步法;Adams外插公式;Adams内插公式;Adams预估—校正公式
1、模型背景及建模
传染病经常在世界各地流行,如霍乱、天花、艾滋病、SARS、H5N1等.建立传染病的数学模型,分析其变化规律,防止其蔓延是一项重要而艰巨的任务.此处仅就一般的传染规律讨论传染病的数学模型.
假设传染病传播期间该地区总人口不变,为常数n.开始染病人数为
,在时刻t的健康人数为
,染病人数为
.由于总人数为n,所以有:
.(1.1)
设单位时间内一个病人能传染的人数与当时健康人数呈正比,比例为常数k,称k为传染系数,于是:
(1.2)
注意到(1.1)式可得:
(1.3)
上述模型就是SI模型,即易感染者(Susceptible)和已感染者(Infective)模型.
2、Adams方法的介绍
Adams方法本质上就是由数值积分方法构造的线性多步法的一种形式.对于非线性常微分方程:
,(2.1)
将
方程两端从
到
积分得:
(2.2)
对
取等距插值节点
,对应的函数值为
.如果
取不同的值,以及
取不同的插值多项式近似,由上式就可以推导出不同的线性多步公式.
1、Adams外插
式子(2.2)中,如果k=0,并选择
作为茶直接点,对函数
作三次差值多项式:
,(2.3)
差值余项为:
(2.4)
把
代入(2.2)得到:
(2.5)
略去第三项,得:
.
对于上式积分部分用变量代换
,并注意:
则:
由以上分析可得线性四步Adams显示公式:
(2.6)
其局部截断误差就是数值积分的误差,即:
(2.7)
因为
在
上不变号,并假设
在
上连续,利用积分中值定理,存在
使得:
(2.7)
因为插值多项式
是在
上作出的,而积分区间为
,故又称为Adams外插公式.
2、Adams内插
若取
,并选择
作为插值节点,作函数
的三次插值多项式.类似于外插公式,可得:
(2.8)
(2.9)
该公式称为Adams内插公式,又称三步隐式方法.
3、Adams预估—校正
由于Adams内插公式是隐式方法,故用它做计算需使用迭代法.通常把Adams外插公式与内插公式结合起来使用,先由前者提供初值,再由后者进行修正,即
(2.10)
当
在求解区域内成立时,迭代收敛.
但若上式中的第二式只迭代一次便得到Adams预估-校正格式
(2.11)
3、对SI模型进行Adams方法求解
对于(1.3)中,k取值为1.1,n取值为13600,
取值为300,区间是(0,1).以下仅使用外插法和预估校正法进行处理.
1、使用Adams外插法
function[z]=f(x,y)
k=1.1;
n=13600;
z=k*y*(n-y)
function[z]=fy(x);
z=13600/(1+(13600/300-1)*exp(-1.1*13600*x));
function[y]=compute(n)
%n是区间剖分的段数;
%x在[a,b]上
a=0;b=1;
h=(b-a)/n;
x=a:
h:
b;
y=zeros(n+1,1);
%x(i)=a+i*h;
%x(i-1)=a+i*h-h;
%x(i-2)=a+i*h-2*h;
%x(i-3)=a+i*h-3*h
y
(1)=0;
y
(2)=fy(h);
y(3)=fy(2*h);
y(4)=fy(3*h);
fori=4:
n
y(i+1)=y(i)+h*(1/24)*(55*f(x(i),y(i))-59*f(x(i-1),y(i-1))+37*f(x(i-2),y(i-2))-9*f(x(i-3),y(i-3)));
end
function[h]=err(n)
a=0;b=1;
h=(b-a)/n;
x=a:
h:
b;
y=compute(n);
fori=1:
(n+1)
realy(i)=fy(x(i));
end
%size(y)
%size(realy)
p=y'-realy
h=sum(abs(y'-realy));
subplot(1,1,1)
plot(x,y,'.',x,realy,'r')
title('显示格式')
运行err函数得
2、使用Adams预估—校正
function[z]=f(x,y)
k=1.1;
n=13600;
z=k*y*(n-y)
function[z]=fy(x);
z=13600/(1+(13600/300-1)*exp(-1.1*13600*x));
function[y]=computew(n)
%n是区间剖分的段数;
%x在[a,b]上
a=0;b=1;
h=(b-a)/n;
x=a:
h:
b;
y=zeros(n+1,1);
y
(1)=0;
y
(2)=fy(h);
y(3)=fy(2*h);
y(4)=fy(3*h);
fori=4:
n
t=y(i)+h*(1/24)*(55*f(x(i),y(i))-59*f(x(i-1),y(i-1))+37*f(x(i-2),y(i-2))-9*f(x(i-3),y(i-3)));
%y(i+1)=(1/(1-9*h/24))*(y(i)+h*(1/24)*(19*f(x(i),y(i))-5*f(x(i-1),y(i-1))+f(x(i-2),y(i-2)))+9*h*exp(x(i+1))/24);
y(i+1)=y(i)+h*(1/24)*(9*f(x(i+1),t)+19*f(x(i),y(i))-5*f(x(i-1),y(i-1))+f(x(i-2),y(i-2)));
end
end
function[h]=errw(n)
a=0;b=1;
h=(b-a)/n;
x=a:
h:
b;
[y]=computew(n);
fori=1:
(n+1)
realy(i)=fy(x(i));
end
%size(y)
%size(realy)
p=y'-realy
h=sum(abs(y'-realy));
subplot(1,1,1)
plot(x,y,'.',x,realy,'r')
title('预估格式')
运行errw函数得
1
2
3
4
5
6
7
8
9
10
11
总
显式格式1e4
-0.0300
-1.2356
-0.9383
-04524
-0.2267
-00818
-0,0663
-0.2964
-1.2636
-0.2397
-0.5281
5.1394e4
预估校正1e-4
-300
0
0
0
0
0
0
0
0
0
0
300
参考文献:
[1]王高雄,周之铭,朱思铭,王寿松.常微分方程[M],第三版.高等教育出版社.2006.
[2]黄云清等.数值计算方法[M],第一版.科学出版社.2009.
[3]李荣华,刘播.微分方程数值解法[M],第四版.高等教育出版社.2009.
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 104753150620 姚鸿泰 关于SI模型的线性多步法求解讨论 关于 SI 模型 线性 步法 求解 讨论