GARCH模型在Matlab中的实现

GARCH模型在Matlab中的实现
GARCH模型在Matlab中的实现

多元GARCH模型预测的Matlab程序

function [parameters, loglikelihood, Ht, likelihoods, stdresid, stderrors, A, B, scores] = full_bekk_mvgarch(data,p,q, BEKKoptions);

% PURPOSE:

% To Estimate a full BEKK multivariate GARCH model. %

%

% USAGE:

% [parameters, loglikelihood, Ht, likelihoods, stdresid, stderrors, A, B, scores] = full_bekk_mvgarch(data,p,q,options);

%

%

% INPUTS:

% data - A t by k matrix of zero mean residuals

% p - The lag length of the innovation process

% q - The lag length of the AR process

% options - (optional) Options for the optimization(fminunc)

%

% OUTPUTS:

% parameters - A (k*(k+1))/2+p*k^2+q*k^2 vector of estimated parameteters. F

% or any k^2 set of Innovation or AR parameters X,

% reshape(X,k,k) will give the correct matrix

% To recover C, use ivech(parmaeters(1:(k*(k+1))/2)

% loglikelihood - The loglikelihood of the function at the optimum

% Ht - A k x k x t 3 dimension matrix of conditional covariances % likelihoods - A t by 1 vector of individual likelihoods

% stdresid - A t by k matrix of multivariate standardized residuals

% stderrors - A numParams^2 square matrix of robust Standad Errors(A^(-1)*B*A^(-1)*t^(-1))

% A - The estimated inverse of the non-robust Standard errors % B - The estimated covariance of teh scores

% scores - A t by numParams matrix of individual scores

% need to try and get some smart startgin values

if size(data,2) > size(data,1)

data=data';

end

[t k]=size(data);

k2=k*(k+1)/2;

scalaropt=optimset('fminunc');

scalaropt=optimset(scalaropt,'TolFun',1e-1,'Display','iter','Diagnostics','on','DiffMaxC hange',1e-2);

startingparameters=scalar_bekk_mvgarch(data,p,q,scalaropt);

CChol=startingparameters(1:(k*(k+1))/2);

C=ivech(startingparameters(1:(k*(k+1))/2))*ivech(startingparameters(1:(k*(k+1))/2))' ;

newA=[];

newB=[];

for i=1:p

newA=[newA diag(ones(k,1))*startingparameters(((k*(k+1))/2)+i)];

end

for i=1:q

newB=[newB diag(ones(k,1))*startingparameters(((k*(k+1))/2)+i+p)];

end

newA=reshape(newA,k*k*p,1);

newB=reshape(newB,k*k*q,1);

startingparameters=[CChol;newA;newB];

if nargin<=3 | isempty(BEKKoptions)

options=optimset('fminunc');

options.Display='iter';

options.Diagnostics='on';

options.TolX=1e-4;

options.TolFun=1e-4;

options.MaxFunEvals=5000*length(startingparameters);

options.MaxIter=5000*length(startingparameters);

else

options=BEKKoptions;

end

parameters=fminunc('full_bekk_mvgarch_likelihood',startingparameters,options,data, p,q,k,k2,t);

[loglikelihood,likelihoods,Ht]=full_bekk_mvgarch_likelihood(parameters,data,p,q,k,k 2,t);

loglikelihood=-loglikelihood;

likelihoods=-likelihoods;

% Standardized residuals

stdresid=zeros(size(data));

for i=1:t

stdresid(i,:)=data(i,:)*Ht(:,:,i)^(-0.5);

end

%Std Errors

if nargout>=6

A=hessian_2sided('full_bekk_mvgarch_likelihood',parameters,data,p,q,k,k2,t);

h=max(abs(parameters/2),1e-2)*eps^(1/3);

hplus=parameters+h;

hminus=parameters-h;

likelihoodsplus=zeros(t,length(parameters));

likelihoodsminus=zeros(t,length(parameters));

for i=1:length(parameters)

hparameters=parameters;

hparameters(i)=hplus(i);

[HOLDER, indivlike] = full_bekk_mvgarch_likelihood(hparameters,data,p,q,k,k2,t);

likelihoodsplus(:,i)=indivlike;

end

for i=1:length(parameters)

hparameters=parameters;

hparameters(i)=hminus(i);

[HOLDER, indivlike] = full_bekk_mvgarch_likelihood(hparameters,data,p,q,k,k2,t);

likelihoodsminus(:,i)=indivlike;

end

scores=(likelihoodsplus-likelihoodsminus)./(2*repmat(h',t,1));

B=cov(scores);

A=A/t;

stderrors=A^(-1)*B*A^(-1)*t^(-1);

end

交通流中的nasch模型及matlab代码元胞自动机

元胞自动机NaSch模型及其MATLAB代码 作业要求 根据前面的介绍,对NaSch模型编程并进行数值模拟: ●模型参数取值:Lroad=1000,p=0.3,Vmax=5。 ●边界条件:周期性边界。 ●数据统计:扔掉前50000个时间步,对后50000个时间步进行统计,需给出的 结果。 ●基本图(流量-密度关系):需整个密度范围内的。 ●时空图(横坐标为空间,纵坐标为时间,密度和文献中时空图保持一致, 画 500个时间步即可)。 ●指出NaSch模型的创新之处,找出NaSch模型的不足,并给出自己的改进思 路。 ●流量计算方法: 密度=车辆数/路长; 流量flux=density×V_ave。 在道路的某处设置虚拟探测计算统计时间T内通过的车辆数N; 流量flux=N/T。 ●在计算过程中可都使用无量纲的变量。 1、NaSch模型的介绍 作为对184号规则的推广,Nagel和Schreckberg在1992年提出了一个模拟车辆交通的元胞自动机模型,即NaSch模型(也有人称它为NaSch模型)。 ●时间、空间和车辆速度都被整数离散化。

● 道路被划分为等距离的离散的格子,即元胞。 ● 每个元胞或者是空的,或者被一辆车所占据。 ● 车辆的速度可以在(0~Vmax )之间取值。 2、NaSch 模型运行规则 在时刻t 到时刻t+1的过程中按照下面的规则进行更新: (1)加速:),1min(max v v v n n +→ 规则(1)反映了司机倾向于以尽可能大的速度行驶的特点。 (2)减速:),min(n n n d v v → 规则(2)确保车辆不会与前车发生碰撞。 (3)随机慢化: 以随机概率p 进行慢化,令:)0, 1-min(n n v v → 规则(3)引入随机慢化来体现驾驶员的行为差异,这样既可以反映随机加速行为, 又可以反映减速过程中的过度反应行为。这一规则也是堵塞自发产生的至关重要因素。 (4)位置更新:n n n v x v +→ ,车辆按照更新后的速度向前运动。 其中n v ,n x 分别表示第n 辆车位置和速度;l (l ≥1)为车辆长度; 11--=+n n n x x d 表示n 车和前车n+1之间空的元胞数;p 表示随机慢化概率;max v 为最大速度。 3、NaSch 模型实例 根据题目要求,模型参数取值:L=1000,p=0.3,Vmax=5,用matlab 软件进行编程,扔掉前11000个时间步,统计了之后500个时间步数据,得到如下基本图和时空图。 3.1程序简介 初始化:在路段上,随机分配200个车辆,且随机速度为1-5之间。 图3.1.1是程序的运行图,图3.1.2中,白色表示有车,黑色是元胞。

Matlab实验报告五(微分方程求解Euler折线法)-推荐下载

数学与信息科学系实验报告 实验名称微分方程求解 所属课程数学软件与实验 实验类型综合型实验 专业信息与计算科学 班级 学号姓名指导教师 吊 顶 到 位 。 连 接 管 半 径 标 式 , 为备 , 查 所 有 试 卷 设 备 进 电 保 护 试 卷 ; 对 整 置 技 术 限 度 内 卷 破 避 免 错 中 资 料

一、实验概述 【实验目的】 熟悉在Matlab 环境下求解常微分方程组和偏微分方程组的方法,掌握利用Matlab 软件进行常微分方程组和偏微分方程组的求解。 【实验原理】 1.dsolve(‘equ1’,’equ2’,...):matlab 求微分方程的解析解。 2.simplify(s):对表达式S 使用MAPLE 的化简规则进行化简。 3.[x,y]=dslove(‘方程1’,‘方程2’,...‘初始条件1’‘初始条件2’,..’自变量’):用字符串方程表示,自变量缺省值为t. 4.ezplot(x,y,[tmin,tmax]):符号函数的作图命令。【实验环境】 MatlabR2010b 二、实验内容 问题1. 求微分方程组在初始条件下的解,并 00dx x y dt dy x y dt ?++=????+-=??00|1,|0t t x y ====[0,0.5]t ∈画出函数的图像. ()y f x =1.分析问题 本题是根据初始条件求微分方程组的特解,并根据t 的范围画出函数的图形。 2.问题求解 syms x y t [x,y]=dsolve('Dx+x+y=0','Dy+x-y=0','x(0)=1','y(0)=0','t')x=simple(x)y=simple(y) ezplot(x,y,[0,0.5]);axis auto 3.结果 x = exp(2^(1/2)*t)/2 + 1/(2*exp(2^(1/2)*t)) - (2^(1/2)*exp(2^(1/2)*t))/4 + 2^(1/2)/(4*exp(2^(1/2)*t)) y = 2^(1/2)/(4*exp(2^(1/2)*t)) - (2^(1/2)*exp(2^(1/2)*t))/4 x = cosh(2^(1/2)*t) - (2^(1/2)*sinh(2^(1/2)*t))/2 、管路敷设技术通过管线不仅可以解决吊顶层配置不规范高中资料试卷问题,而且可保障各类管路习题到位。在管路敷设过程中,要加强看护关于管路高中资料试卷连接管口处理高中资料试卷弯扁度固定盒位置保护层防腐跨接地线弯曲半径标高等,要求技术交底。管线敷设技术包含线槽、管架等多项方式,为解决高中语文电气课件中管壁薄、接口不严等问题,合理利用管线敷设技术。线缆敷设原则:在分线盒处,当不同电压回路交叉时,应采用金属隔板进行隔开处理;同一线槽内,强电回路须同时切断习题电源,线缆敷设完毕,要进行检查和检测处理。、电气课件中调试对全部高中资料试卷电气设备,在安装过程中以及安装结束后进行 高中资料试卷调整试验;通电检查所有设备高中资料试卷相互作用与相互关系,根据生产工艺高中资料试卷要求,对电气设备进行空载与带负荷下高中资料试卷调控试验;对设备进行调整使其在正常工况下与过度工作下都可以正常工作;对于继电保护进行整核对定值,审核与校对图纸,编写复杂设备与装置高中资料试卷调试方案,编写重要设备高中资料试卷试验方案以及系统启动方案;对整套启动过程中高中资料试卷电气设备进行调试工作并且进行过关运行高中资料试卷技术指导。对于调试过程中高中资料试卷技术问题,作为调试人员,需要在事前掌握图纸资料、设备制造厂家出具高中资料试卷试验报告与相关技术资料,并且了解现场设备高中资料试卷布置情况与有关高中资料试卷电气系统接线等情况,然后根据规范与规程规定,制定设备调试高中资料试卷方案。 、电气设备调试高中资料试卷技术电力保护装置调试技术,电力保护高中资料试卷配置技术是指机组在进行继电保护高中资料试卷总体配置时,需要在最大限度内来确保机组高中资料试卷安全,并且尽可能地缩小故障高中资料试卷破坏范围,或者对某些异常高中资料试卷工况进行自动处理,尤其要避免错误高中资料试卷保护装置动作,并且拒绝动作,来避免不必要高中资料试卷突然停机。因此,电力高中资料试卷保护装置调试技术,要求电力保护装置做到准确灵活。对于差动保护装置高中资料试卷调试技术是指发电机一变压器组在发生内部故障时,需要进行外部电源高中资料试卷切除从而采用高中资料试卷主要保护装置。

欧拉法matlab程序

法 function [x,y]=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end x=x';y=y'; x1=0::1;y1=(1+2*x1).^; plot(x,y,x1,y1) >> dyfun=inline('y-2*x/y'); [x,y]=naeuler(dyfun,[0,1],1,;[x,y] ans = 2.隐式Euler法 function [x,y]=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h); end x=x';y=y'; x1=0::1;y1=(1+2*x1).^; plot(x,y,x1,y1) function y=iter(dyfun,x,y,h) y0=y;e=1e-4;K=1e+4; y=y+h*feval(dyfun,x,y); y1=y+2*e;k=1; while abs(y-y1)>e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1; if k>K error('迭代发散'); end end >> dyfun=inline('y-2*x/y');

[x,y]=naeulerb(dyfun,[0,1],1,;[x,y] ans = 3.改进Euler法 function [x,y]=naeuler2(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*k1; k2=feval(dyfun,x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2; end x=x';y=y'; x1=0::1;y1=(1+2*x1).^; plot(x,y,x1,y1) >> dyfun=inline('y-2*x/y'); [x,y]=naeuler2(dyfun,[0,1],1,;[x,y] ans =

实验一 用MATLAB处理系统数学模型

实验一用MATLAB处理系统数学模型 一、实验原理 表述线性定常系统的数学模型主要有微分方程、传递函数、动态结构图等.求拉氏变换可用函数laplace(ft,t,s),求拉式反变换可用函数illaplace(Fs,s,t);有关多项式计算的函数主要有roots(p),ploy(r),conv(p,q),ployval(n,s);求解微分方程可采用指令 s=dslove(‘a_1’,’a_2’,’···,’a_n’);建立传递函数时,将传递函数的分子、分母多项式的系数写成两个向量,然后用tf()函数来给出,还可以建立零、极点形式的传递函数,采用的函数为zpk(z,p,k);可用函数sys=series(sys1,sys2)来实现串联,用 sys=parallel(sys1,sys2)来实现并联,可用函数sys=feedback(sys1,sys2,sign)来实现系统的反馈连接,其中sign用来定义反馈形式,如果为正反馈,则sign=+1,如果为负反馈,则sign=-1。 二、实验目的 通过MATLAB软件对微分方程、传递函数和动态结构图等进行处理,观察并分析实验结果。 三、实验环境 MATLAB2012b 四、实验步骤 1、拉氏变换 syms s t; ft=t^2+2*t+2; st=laplace(ft,t,s) 2、拉式反变换 syms s t; Fs=(s+6)/(s^2+4*s+3)/(s+2); ft=ilaplace(Fs,s,t) 3、多项式求根 p=[1 3 0 4]; r=roots(p) p=poly(r) 4、多项式相乘 p=[ 3 2 1 ];q=[ 1 4];

欧拉图fluery算法matlab

clear all A=zeros(16); for i=1:16 %A(i,i)=1/2; if i+1<=16 && mod(i,4)~=0 A(i,i+1)=1; end if i+4<=16 A(i,i+4)=1; end end A(2,5)=1; A(3,8)=1; A(9,14)=1; A(12,15)=1; % A(1,6)=1; % A(6,11)=1; % A(11,16)=1; % A(16,1)=1; A=A+A'; [T3 c3]=Fleuf2(A); pos3(1:2,1)=0; for i=1:16 if i==1, pos3(1:2,i)=0; else if mod(i-1,4)~=0 pos3(1,i)=pos3(1,i-1)+1; pos3(2,i)=pos3(2,i-1); else pos3(1,i)=pos3(1,i-4); pos3(2,i)=pos3(2,i-4)-1; end end end figure (1), title('3rd Question')

for j=i:16 if A(i,j)==1, plot([pos3(1,i),pos3(1,j)],[pos3(2,i),pos3(2,j)]); hold on, end end end for i=2:c3 draw_arrow(pos3(:,T3(1,i))',pos3(:,T3(2,i))',0.5) x = mean(pos3(1,T3(:,i))); y = mean(pos3(2,T3(:,i))) text(x,y,num2str(i),'FontSize',18); pause; end pause off; hold off function [T c]=Fleuf1(d) n = length(d); b = d; b(b == Inf)=0; b(b~=0) = 1; m = 0; a = sum(b); eds = sum(a)/2; ed = zeros(2,eds); vexs = zeros(1,eds+1); matr = b; for i=1:n if mod(a(i),2) ==1 m = m+1; end end if m~=0 fprintf('there is not exist Euler path.\n'); T=0;c=0; end if m==0

用matlab实现碰撞模型程序代码

用m a t l a b实现碰撞模型程序代码 标准化工作室编码[XX968T-XX89628-XJ668-XT689N]

c l c; clear; fill([6,7,7,6],[5,5,0,0],[0,0.5,0]);%右边竖条的填充 holdon;%保持当前图形及轴系的所有特性 fill([2,6,6,2],[3,3,0,0],[0,0.5,0]);%左边竖条的填充 holdon;%保持当前图形及轴系的所有特性 t1=0:pi/60:pi; plot(4-2*sin(t1-pi/2),5-2*cos(t1-pi/2));%绘制中间的凹弧图形gridon;%添加网格线 axis([0,9,0,9]);%定义坐标轴的比例% axis('off');%关闭所有轴标注,标记,背景 fill([1,2,2,1],[5,5,0,0],[0,0.5,0]);%中间长方形的填充 holdon;%保持当前图形及轴系的所有特性 title('碰撞');%定义图题 x0=6; y0=5; head1=line(x0,y0,'color','r','linestyle','.','erasemode','xor','marke rsize',30); head2=line(x0,y0,'color','r','linestyle','.','erasemode','xor','marke rsize',50);%设置小球颜色,大小,线条的擦拭方式 t=0;%设置小球的初始值 dt=0.001;%设置运动周期 t1=0;%设置大球的初始值 dt1=0.001; while1%条件表达式 t=t+dt; x1=9-1*t; y1=5; x3=6; y3=5; ift>0 x2=6; y2=5;%设置小球的运动轨迹 end ift>2.8 t=t+dt; a=sin(t-3); x1=6.1; y1=5.1; x3=4-2*sin(1.5*a); y3=5-2*cos(1.5*a);%设置大球的运动轨迹 end

欧拉法matlab程序学习课件.doc

1.Euler法 function [x,y]=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end x=x';y=y'; x1=0:0.2:1;y1=(1+2*x1).^0.5; plot(x,y,x1,y1) >> dyfun=inline('y-2*x/y'); [x,y]=naeuler(dyfun,[0,1],1,0.2);[x,y] ans = 0 1.0000 0.2000 1.2000 0.4000 1.3733 0.6000 1.5315 0.8000 1.6811 1.0000 1.8269 2.隐式Euler法 function [x,y]=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h); end x=x';y=y'; x1=0:0.2:1;y1=(1+2*x1).^0.5;

plot(x,y,x1,y1) function y=iter(dyfun,x,y,h) y0=y;e=1e-4;K=1e+4; y=y+h*feval(dyfun,x,y); y1=y+2*e;k=1; while abs(y-y1)>e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1; if k>K error('迭代发散'); end end >> dyfun=inline('y-2*x/y'); [x,y]=naeulerb(dyfun,[0,1],1,0.2);[x,y] ans = 0 1.0000 0.2000 1.1641 0.4000 1.3014 0.6000 1.4146 0.8000 1.5019 1.0000 1.5561 3.改进Euler法 function [x,y]=naeuler2(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n)); y(n+1)=y(n)+h*k1; k2=feval(dyfun,x(n+1),y(n+1));

飞机碰撞模型

飞机碰撞模型 摘要 第六架在边长为160km的正方形区域内以的飞行角从坐标为(0,0)的点出发,在飞行过程中不与其它五架飞机发生碰撞,即在该区域内与其它任意飞机的距离大于8km,就要不断调整该飞机的飞行角度,使其任意时刻与其他飞机的距离大于8km,利用空间中点的距离定义,计算任意时刻该飞机与其他飞机的距离,找到调整角度的最小值为。 1、问题重述 在约10000km高空的某边长160km的正方形区域内,有5架飞机均以800km/h的速度作水平飞行,不碰撞的标准为在该区域内任意两架飞机的距离大于8km。现有5架飞机在区域内飞行且它们不会碰撞,其初始坐标和飞行方向由下表给出: 现有第6架飞机要进入该区域,坐标为(0,0),飞行角为,如果其与内部的5架飞机发生碰撞,就需要调整其飞行角度,请建立优化模型,确定其与内部5架飞机不碰撞的最小调整角。 2、基本假设 1、五架飞机在规定正方形区域飞行中不随意改变路线; 2、飞机在飞行中不考虑其他未知因素; 3、符号说明 :正方形区域的边长; :第i架飞机飞行的方向角度; :第六架飞机飞行过程中的调整角度; :第架、第架飞机的距离; :第架飞机在区域内飞行的路线长度; :第架飞机的飞行速度; :第架飞机在区域内的飞行时间; :第i架飞机的横坐标; :第i架飞机的纵坐标; 4、模型的建立与求解 1、模型的建立 先根据五架飞机起始点与终点坐标,在规定的网格区域内画出它们的飞行路线,再根据给出的区域长度与各架飞机飞行速度,计算出各架飞机在区域内的飞行时间, 再根据计算得出的时间,得出时刻各架飞机的坐标,求出在该时刻第六架飞机与其他五架飞机的距离 即 当<8时,此时就需要调整第六架飞机的飞行角度,使其与另外五架飞机

用matlab实现碰撞模型程序代码

clc; clear; fill([6,7,7,6],[5,5,0,0],[0,0.5,0]);%右边竖条的填充 hold on; %保持当前图形及轴系的所有特性 fill([2,6,6,2],[3,3,0,0],[0,0.5,0]);%左边竖条的填充 hold on;% 保持当前图形及轴系的所有特性 t1=0:pi/60:pi; plot(4-2*sin(t1-pi/2),5-2*cos(t1-pi/2));%绘制中间的凹弧图形 grid on;%添加网格线 axis([0,9,0,9]);%定义坐标轴的比例% axis('off');%关闭所有轴标注,标记,背景 fill([1,2,2,1],[5,5,0,0],[0,0.5,0]);%中间长方形的填充 hold on;% 保持当前图形及轴系的所有特性 title('碰撞');%定义图题 x0=6; y0=5; head1=line(x0,y0,'color','r','linestyle','.','erasemode','xor','markersize',30); head2=line(x0,y0,'color','r','linestyle','.','erasemode','xor','markersize',50); %设置小球颜色,大小,线条的擦拭方式 t=0;%设置小球的初始值 dt=0.001;%设置运动周期 t1=0;%设置大球的初始值 dt1=0.001; while 1%条件表达式 t=t+dt; x1=9-1*t; y1=5; x3=6; y3=5; if t>0 x2=6; y2=5;%设置小球的运动轨迹 end if t>2.8 t=t+dt; a=sin(t-3); x1=6.1; y1=5.1; x3=4-2*sin(1.5*a); y3=5-2*cos(1.5*a);%设置大球的运动轨迹

实验五 欧拉法Matlab实验报告

北京理工大学珠海学院实验报告 ZHUHAI CAMPAUS OF BEIJING INSTITUTE OF TECHNOLOGY 班级2012电气2班学号120109021010姓名陈冲指导教师张凯成绩 实验题目(实验五)欧拉法实验地点及时间JD501 2014/1/2(6-7节) 一、实验目的 1.掌握用程序语言来编辑函数。 2.学会用MATLAB编写Euler.m以及TranEuler.m函数。 二、实验环境 Matlab软件 三、实验内容 1、以书中第124页题目11为例编辑程序来实现计算结果。 2、使用MATLAB进行编写: 第一步:编写Euler.m函数,代码如下 编写TranEuler.m函数,代码如下 第二步:利用上述函数编辑命令:(可见实验结果中的截图)

在此之前先建立一个名为f.m 的M 文件,代码如下 function z=f(x); z=8-3y; 再编辑代码: 得到了欧拉法的结果:y (0.4)=2.47838030901267 编辑另一段命令: 得到改进欧拉法的结果:y (0.4)=2.46543714659780 在此基础上,我还编辑龙格库达的命令窗口代码,如下: 四、实验题目 用欧拉法和改进欧拉法求解初值问题'83,(0)2y y y =-=,试取步长0.2h =计算(0.4)y 的近似值。 五、实验结果

六、总结 通过这次实验我掌握了将得到的解进一步精确,而且要学会比较这几种方法的精确性,显然,四阶龙格库达比改进欧拉发精确,改进欧拉发比欧拉法精确。 实验难度不大,要比较n的取值不同,产生的影响不同。

基于MATLAB的地震正演模型实现[1]

基于MATLAB的地震正演模型实现 贾跃玮 (中国地质大学(北京) 北京100083) 摘 要 人工合成地震正演模型是进行三维模型计算的基础。针对地震勘探的原理,本文运用MATLAB强大数学计算和图像可视化功能,对一个三层介质模型制作了人工合成地震记录。文章首先说明了地震记录形成的物理机制,然后介绍了地质模型的构造及参数选择,最后针对该具体地质模型制作了合成地震记录。 关键词 地震;MATLAB;正演 0引 言 地震勘探就是利用地下介质弹性和密度的差异,通过观测和分析大地对人工激发地震波的响应,推断地下岩层的性质和形态的地球物理方法。地震勘探是钻探前勘测石油与天然气资源的重要手段,在煤田和工程地质勘查、区域地质研究和地壳研究等方面,也得到广泛应用。 人工合成二维地震模型记录是各种复杂地震模型正演计算的基础,是对地震勘探经典理论的忠实实现。在实际工作中,针对具体地质构造进行二维地震模拟能够有效帮助地球物理工作者在地震剖面上识别各种地质现象。MATLAB环境集编程、画图于一体,特别适合人工合成地震记录的快速实现。因此,我们在MATLAB环境下设计了一个三层地质模型,并对该模型模拟了地震记录,旨在可视化地观察地震波场记录特征并验证地震褶积模型。 1地震记录形成的物理机制 在地震记录上看到的波形是地震子波叠加的结果,从地下许多反射界面发生反射时形成的地震子波,振幅大小决定于反射界面反射系数的绝对值,极性的正负决定于反射系数的正负,到达时间的先后取决于界面深度和覆盖层的波速。若地震子波波形用S(t)表示,反射系数是双程垂直反射旅行时t的函数,用R(t)表示,地震记录f(t)形成的物理过程在数学上就可表示为:f(t)=S(t)3R(t)=∫0T S(τ)R(t-τ)dτ 地震子波和反射系数资料常常不易取得,因此计算时常做这样一些假设: (1)地质模型的建立是来自大量观察实际地质结构的经验性归纳总结。 (2)为了模型建立和计算过程中突出理论数值,去除了一些干扰因素,对一切衰减、噪声都不进行考虑。 (3)地层在横向上均匀,纵向上是由大量具有不同弹性性质的薄层构成。 (4)地震子波以平面波形式垂直入射到界面,各薄层的反射子波与地震子波形状相同,只是振幅及极性不同。 (5)所有波的转换、吸收及绕射等能量损失都不考虑。 基于以上这些假设条件进行地震记录合就必须已知地震子波以及地层的反射系数,而反射系数又主要由地层的波阻抗反映,所以必须首先获取地层的速度和密度资料。 速度资料可通过连续速度测井获得,密度资料可从密度测井获得,得不到密度资料时,可近似假定密度不变,以速度曲线代替波阻抗曲线来计算反射系数。加德纳根据实际资料提出了一个由速度推算密度的经验公式: ρ=0.23V0.25 (速度单位:英尺/秒) 或 ρ=0.31V0.25 (速度单位:米/秒)

matlab 欧拉算法 附截图

设系统方程为:y t y y /2)1(-=,1)0(=y ,用改进欧拉法求解各离散点y 的数值解,步长 10,1.0≤≤=t h ,解析解为t y 21+= 。 解:改进欧拉法 ),(1n n n p n y t hf y y +=+ )],(),([5.0111p n n n n n c n y t f y t f h y y +++++= 已知 n n n n n y t y y t f /2),(-= n n n n n n n p n y ht y h y t y h y y /2)1()/2(1-+=-+=+ 1 111111/5.0/)5.01()]/2()/2[(5.0+++++++-+-+=-+-+=n n n n n n n n n n n n n c n y ht hy y ht y h y t y y t y h y y 程序: h=0.1; t=0:h:1; N=length(t); y=ones(1,N); ey=ones(1,N); zy=ones(1,N); for k=1:N-1 y(1,k+1)=(1+h)*y(1,k)-(2*h*(k-1)/(N-1))./y(1,k);%预估公式 ey(1,k+1)=(1+h)*ey(1,k)-(2*h*(k-1)/(N-1))./ey(1,k);%欧拉公式 y(1,k+1)=(1+0.5*h)*y(1,k)-(h*(k-1)/(N-1))./y(1,k)+0.5*h*y(1,k+1)-(h*k/(N-1))./y(1,k+1);%改进欧拉 zy(1,k+1)=(1+2*k/(N-1)).^0.5;%解析解 end plot(t,zy,'-xk',t,y,':ob',t,ey,'-.*r','linewidth',1.0); xlabel('t'); ylabel('y'); 截图:

Euler法解微分方程-Matlab程序

%主程序main.m-----OK! clear; T=0.0; Y=zeros(3,1); Y(1)=1.0;Y(2)=1.0;Y(3)=1.0; H1=0.05;M=3;EPS=1.0e-05;%EPS精度要求M方程个数H1拟定的输出步长 for i=1:10 [X,Y]=euler1(T,H1,Y,M,EPS) T=T+H1; End %变步长euler方法 function [X,Y1]=euler1(T,H1,Y,M,EPS) %M-方程个数,EPS-精度,Y0-右端初值,T-自变量前一点值,H-步长 N=1;P=1+EPS;X=T;G=zeros(M,1); H=H1;%H-在程序中要改变的步长H1-主程序中确定的输出步长 for i=1:M C(i)=Y(i); end K1=zeros(M,1);K2=zeros(M,1);K3=zeros(M,1);K4=zeros(M,1); while P>=EPS %变步长积分一步(H1) for i=1:M G(i)=Y(i); Y(i)=C(i); end DT=H/N; T=X; %--变步长积分过程 for j=1:N K1=F(Y); K2=F(Y+H/2*K1'); K3=F(Y+H/2*K2'); K4=F(Y+H*K3'); for i=1:M Y(i)=Y(i)+H/6*(K1(i)+2*K2(i)+2*K3(i)+K4(i)); T=T+DT; end end %--------------------- P=0.0; for i=1:M Q=abs(Y(i)-G(i)); if Q>P

P=Q; end end H=H/2.0; N=N+N; end T=X; X=T+H1; Y1=Y; %右端函数值function D=F(y) D(1)=y(2); D(2)=-1*y(1); D(3)=y(3);

完全弹性碰撞matlab

Matlab设计实验 课题名称:完全弹性碰撞 一.设计背景: 完全弹性碰撞(Perfect Elastic Collision):在理想情况下,完全弹性碰撞的物理过程满足动量守恒和能量守恒。如果两个碰撞小球的质量相等,联立动量守恒和能量守恒方程时可解得:两个小球碰撞后交换速度。如果被碰撞的小球原来静止,则碰撞后该小球具有了与碰撞小球一样大小的速度,而碰撞小球则停止。多个小球碰撞时可以进行类似的分析。 二.设计意义 真实情况下,由于小球间的碰撞并非理想的弹性碰撞,还会有能量的损失,所以最后小球还是要停下来。 所以该设计主要用于研究能量守恒中的某些问题。还有就是用于实验演示。三.程序设计 该程序主要设置了三个不同颜色的小球,在真空环境下(理想环境下)的碰撞实验演示。 该程序可以通过改变各种参数,研究各种情况下的实验数据。 程序: pole=1.8;%定义摆线的长度 xmax=2;%定义横坐标长度 ymax=2;%定义纵坐标长度 basew=2.3;%定义图中方框的宽度 baseh=2.3;%定义图中方框的高度 instant=0.2;%定义摆线间距 %三视图的初始设置 %第一幅图

figure('name','理想情况下能量守恒定律 1','position',[500,340,440,320]);%定义第一幅图的标题和位置 fill([xmax,xmax,-xmax,-xmax,xmax,xmax-0.05,xmax-0.05,- xmax+0.05,-xmax+0.05,xmax-0.05],[ymax,-ymax,- ymax,ymax,ymax,ymax-0.05,-ymax+0.05,-ymax+0.05,ymax- 0.05,ymax-0.05],[0,1,1]); %填充底座背景 hold on;%保持当前图形及坐标所有特性 fill([xmax-0.05,xmax-0.05,-xmax+0.05,-xmax+0.05],[ymax- 0.5 ,ymax-0.55,ymax-0.55,ymax-0.5],'g');%填充方框内横杆背景 hold on;%保持当前图形及坐标所有特性 text(-0.25,1.7,'1');text(0,1.7,'2');text(0.25,1.7,'3');%在坐标处标识 说明文字 text( -1.0,1.7,'a');text( -1.0,-1.7,'b');%在坐标处标识说明文字 text(1.0,1.7,'真空容器');text(-1.8,1.7,'主视图');%在坐标处标识说明文 字 axis([-basew,basew,-baseh,baseh]);%定义背景坐标范围在x(-2.3~2.3) Y(-2.3~2.3)之间 %axis('off');%覆盖坐标刻度并填充背景 theta0=7 *pi/6;%摆线1的初始角度 x0=pole*cos(theta0);%摆线1末端x坐标 y0=pole*sin(theta0)+1.5;%摆线1末端y坐标 body1=line([-instant,x0-instant],[1.5,y0],'color','r','linestyle','- ','erasemode','xor');%设置摆线1 head1=line(x0- instant,y0,'color','r','linestyle','.','erasemode','xor','markersize',40);%设置第一个小球颜色,大小 theta1=3*pi/2;%摆线2,3的角度 x1=pole*cos(theta1);%摆线2,3末端x坐标 y1=pole*sin(theta1)+1.5;%摆线2,3末端y坐标 body=line([-0.001,x1],[1.5,y1],'color','k','linestyle','- ','erasemode','xor');%设置摆线2

MATLABEuler法解常微分方程

Euler法解常微分方程 Euler法解常微分方程算法: Step 1 分别取积分上限、积分下限、步长 Step 2计算判断是否成立,成立转到Step 3,否则继续进行Step 4 Step 3 计算 Step 4 Euler法解常微分方程算程序: function euler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 x=min(A); b=max(A); y=y0; while x

Step 3 (1)做显性Euler预测 (2)将带入 Step 4计算判断是否成立,成立返回Step 3,否则继续进行Step 5 Step 5 改进Euler法解常微分方程算程序: function gaijineuler2(fun,y0,A,h) %fun--y' %y0---初值 %A----x取值范围 %a----x左区间端点值 %b----x右区间端点值 %h----给定步长 a=min(A); b=max(A); x=a:h:b; y(1)=y0; for i=1:length(x)-1 w1=feval(fun,x(i),y(i)); y(i+1)=y(i)+h*w1; w2=feval(fun,x(i+1),y(i+1)); y(i+1)=y(i)+h*(w1+w2)/2; end x=x' y=y' 例:用改进Euler法计算下列初值问题(取步长h=0.25) 输入:fun=inline('-x*y^2') gaijineuler2(fun,2,[0 5],0.25) 得到: x = 0.2500 0.5000 0.7500 1.0000 1.2500 1.5000 1.7500 2.0000 2.2500 2.5000 2.7500

实验四 用MATLAB求解状态空间模型

实验四 用MATLAB 求解状态空间模型 1、实验设备 MATLAB 软件 2、实验目的 ① 学习线性定常连续系统的状态空间模型求解、掌握MATLAB 中关于求解该模型的主要函数; ② 通过编程、上机调试,进行求解。 3、实验原理说明 Matlab 提供了非常丰富的线性定常连续系统的状态空间模型求解(即系统运动轨迹的计算)的功能,主要的函数有: 初始状态响应函数initial()、阶跃响应函数step()以及可计算任意输入的系统响应数值计算函数lsim()和符号计算函数sym_lsim()。 数值计算问题可由基本的Matlab 函数完成,符号计算问题则需要用到Matlab 的符号工具箱。 4、实验步骤 ① 根据所给状态空间模型,依据线性定常连续系统状态方程的解理论,采用MATLAB 编程。 ② 在MATLAB 界面下调试程序,并检查是否运行正确。 习题1:试在Matlab 中计算如下系统在[0,5s]的初始状态响应,并求解初始状态响应表达式。 Matlab 程序如下: A=[0 1; -2 -3]; B=[]; C=[]; D=[]; x0=[1; 2]; sys=ss(A,B,C,D); [y,t,x]=initial(sys,x0,0:5); plot(t,x) 0011232????==????--????x x x

习题2:试在Matlab 中计算如下系统在[0,10s]内周期为3s 的单位方波输入下的状态响应。并计算该系统的单位阶跃状态响应表达式。 Matlab 程序如下: A=[0 1; -2 -3]; B=[0; 1]; C=[]; D=[]; x0=[1; 2]; sys=ss(A,B,C,D); [u t]=gensig('square',3,10,0.1) 0011232????==????--???? x x x

解微分方程欧拉法-R-K法及其MATLAB实例

解微分方程的欧拉法,龙格-库塔法及其MATLAB简单实例 欧拉方法(Euler method)用以对给定初值的常微分方程(即初值问题)求解? 分为前进EULER法、后退EULER法、改进的EULER法。?缺点:?欧拉法简单地取 切线的端点作为下一步的起点进行计算,当步数增多时,误差会因积累而越来越大。因此欧拉格式一般不用于实际计算。?改进欧拉格式:?为提高精度,需要在 欧拉格式的基础上进行改进。采用区间两端的斜率的平均值作为直线方程的斜率。改进欧拉法的精度为二阶。?算法为: 微分方程的本质特征是方程中含有导数项,数值解法的第一步就是设法消除其导数值。 对于常微分方程: x∈[a,b] y(a) = y0 可以将区间[a,b]分成n段,那么方程在第xi点有y'(xi) = f(xi,y (xi)),再用向前差商近似代替导数则为: 在这里,h是步长,即相邻两个结点间的距离。因此可以根据xi点和yi点的数值计算出yi+1来: i=0,1,2,L 这就是向前欧拉格式。 改进的欧拉公式: 将向前欧拉公式中的导数f(xi,yi)改为微元两端导数的平均,即 上式便是梯形的欧拉公式。

可见,上式是隐式格式,需要迭代求解。为了便于求解,使用改进的欧拉公式: 数值分析中,龙格-库塔法(Runge-Kutta)是用于模拟常微分方程的解的重要的一类隐式或显式迭代法。 实际上,龙格-库塔法是欧拉方法的一种推广,向前欧拉公式将导数项简单取为f(xn,yn),而改进的欧拉公式将导数项取为两端导数的平均。 龙格-库塔方法的基本思想: 在区间[xn,xn+1]内多取几个点,将他们的斜率加权平均,作为导数的近似。 龙格库塔法的家族中的一个成员如此常用,以至于经常被称为“RK4”或者就是“龙格库塔法”。? 令初值问题表述如下。 则,对于该问题的RK4由如下方程给出: 其中 这样,下一个值(yn+1)由现在的值(y n)加上时间间隔(h)和一个估算的斜率的乘积决定。该斜率是以下斜率的加权平均: k1是时间段开始时的斜率;? k2是时间段中点的斜率,通过欧拉法采用斜率k1来决定y在点tn + h/2的值;? k3也是中点的斜率,但是

欧拉法matlab程序资料

程l欧序batam法拉.精品文档 1.Euler法 function [x,y]=naeuler(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n)); end x=x';y=y'; x1=0:0.2:1;y1=(1+2*x1).^0.5;plot(x,y,x1,y1) >> dyfun=inline('y-2*x/y'); [x,y]=naeuler(dyfun,[0,1],1,0.2);[x,y] ans = 0 1.0000 0.2000 1.2000 0.4000 1.3733 0.6000 1.5315

0.8000 1.6811 1.0000 1.8269 2.隐式Euler法 function [x,y]=naeulerb(dyfun,xspan,y0,h) x=xspan(1):h:xspan(2); y(1)=y0; for n=1:length(x)-1 y(n+1)=iter(dyfun,x(n+1),y(n),h); end 收集于网络,如有侵权请联系管理员删除. 精品文档 x=x';y=y'; x1=0:0.2:1;y1=(1+2*x1).^0.5;plot(x,y,x1,y1)function y=iter(dyfun,x,y,h) y0=y;e=1e-4;K=1e+4; y=y+h*feval(dyfun,x,y); y1=y+2*e;k=1; while abs(y-y1)>e y1=y; y=y0+h*feval(dyfun,x,y); k=k+1; if k>K error('迭代发散'); end end >> dyfun=inline('y-2*x/y'); [x,y]=naeulerb(dyfun,[0,1],1,0.2);[x,y] ans = 0 1.0000 0.2000 1.1641 0.4000 1.3014 0.6000 1.4146 0.8000 1.5019 1.0000 1.5561 3.改进Euler法

相关文档
最新文档