简单随机微分方程数值解

news2024/11/15 2:04:55

1.随机微分方程求解:dX(t) =− αXtdt + σdWt

法一:Euler-Maruyama
在这里插入图片描述

%%
%O-U过程
%dX(t)=-alpha*Xt*dt+sigma*dWt,X|t=0=X0
%alpha=2,sigma=1,X0=1
% 设置初始参数
T = 1;                 % 时间区间长度
N = 1000;              % 离散化的时间步数
dt = T/N;              % 时间步长
X = zeros(1,N+1);      % 存储解向量
X(1) = 1;              % 初始条件
alpha = 2;
sigma = 1;

% 模拟数值解
for i = 1:N
    dW = sqrt(dt)*randn();       % 标准正态分布增量
    X(i+1) = X(i) - alpha*X(i)*dt + sigma*dW;   % 欧拉方法更新
end

% 绘图
plot(linspace(0,T,N+1),X,'*-')   % 根据时间步长将x轴离散化
xlabel('t')
ylabel('X(t)')
title('随机微分方程数值解')
legend('Euler数值解')

在这里插入图片描述
法二:Milstein
在这里插入图片描述
在这里插入图片描述

% 设置参数
alpha = 2;
sigma = 1;
X0 = 1;
T = 1;
N = 1000;
dt = T/N;

% 初始化向量和随机项
t = 0:dt:T;
W = [0,cumsum(randn(1,N).*sqrt(dt))];
%使用cumsum函数生成一个与t等长的Wiener过程(随机项)

% 初始化X
X = zeros(1,N+1);
X(1) = X0;

% Milstein方法计算X
for i = 1:N
    dW = W(i+1) - W(i);
    X(i+1) = X(i) - alpha*X(i)*dt + sigma*X(i)*dW ...
             + 0.5*sigma^2*X(i)*(dW^2-dt);
    %在Milstein方法中,我们需要对二次变差数进行估计
    %因此在计算时需要添加正交项0.5*sigma^2*X(i)*(dW^2-dt)
end

% 绘制图形
plot(t,X,'*-')
title('随机微分方程数值解')
xlabel('t')
ylabel('X(t)')
legend('Milstein数值解')

在这里插入图片描述

2.受高斯白噪声激励的系统,FPK求解

考虑一个由下列方程支配的随机激励的单自由度振子:
X ¨ + h ( X , X ˙ ) = X W 1 ( t ) + X ˙ W 2 ( t ) + W 3 ( t ) \ddot{X}+h(X,\dot{X})=XW_{1}(t)+\dot XW_{2}(t)+W_{3}(t) X¨+h(X,X˙)=XW1(t)+X˙W2(t)+W3(t)
其中, h ( X , X ˙ ) h(X,\dot X) h(X,X˙)表示阻尼力和恢复力, W l ( t ) , l = 1 , 2 , 3 W_{l}(t),l=1,2,3 Wl(t),l=1,2,3是独立的高斯白噪声,其相关函数为 E [ W l ( t ) W s ( t + τ ) ] = 2 π K l s δ l s δ ( τ ) E[W_{l}(t)W_{s}(t+\tau)]=2\pi K_{ls} \delta_{ls}\delta(\tau) E[Wl(t)Ws(t+τ)]=2πKlsδlsδ(τ)
这里, δ l s \delta_{ls} δls δ ( τ ) \delta(\tau) δ(τ)不一样,前者是克罗内克(Kronecker) δ \delta δ,即 δ l s = 1 , l = s ; δ l s = 0 , l ≠ s . \delta_{ls}=1,l=s;\delta_{ls}=0,l\neq s. δls=1,l=s;δls=0,l=s. δ ( τ ) \delta(\tau) δ(τ)是狄拉克函数, δ = ∞ , τ = 0 ; δ = 0 , τ ≠ 0. \delta=\infty, \tau=0;\delta=0,\tau \neq 0. δ=,τ=0;δ=0,τ=0.

X 1 = X , X 2 = X ˙ X_{1}=X,X_{2}=\dot{X} X1=X,X2=X˙可转换状态空间的两个方程
X ˙ 1 = X 2 \dot{X}_{1}=X_{2} X˙1=X2
X ˙ 2 = − h ( X 1 , X 2 ) + X 1 W 1 ( t ) + X 2 W 2 ( t ) + W 3 ( t ) \dot{X}_{2}=-h(X_{1},X_{2})+X_{1}W_{1}(t)+X_{2}W_{2}(t)+W_{3}(t) X˙2=h(X1,X2)+X1W1(t)+X2W2(t)+W3(t)

对于 d d t X ( t ) = f ( X , t ) + g ( X , t ) W ( t ) \frac{d}{dt}X(t)=f(X,t)+g(X,t)W(t) dtdX(t)=f(X,t)+g(X,t)W(t)
相应的FPK方程为 ∂ ∂ t p + ∑ j = 1 n ∂ ∂ x j ( a j p ) − 1 2 ∑ j , k = 1 n ∂ 2 ∂ x j ∂ x k ( b j k p ) = 0 \frac{\partial }{\partial t}p+\sum_{j=1}^{n}\frac{\partial }{\partial x _{j}}(a_{j}p)-\frac{1}{2}\sum_{j,k=1}^{n}\frac{\partial ^{2}}{\partial x_{j}\partial x_{k}}(b_{jk}p)=0 tp+j=1nxj(ajp)21j,k=1nxjxk2(bjkp)=0
其中,一、二阶导数矩为 a j ( x , t ) = f j ( x , t ) + π ∑ r = 1 n ∑ l , s = 1 m K l s g r s ( x , t ) ∂ ∂ x r g j l ( x , t ) a_{j}(\mathbf{x},t)=f_{j}(\mathbf{x},t)+\pi\sum_{r=1}^{n}\sum_{l,s=1}^{m}K_{ls}g_{rs}(\mathbf{x},t)\frac{\partial }{\partial x_{r}}g_{jl}(\mathbf{x},t) aj(x,t)=fj(x,t)+πr=1nl,s=1mKlsgrs(x,t)xrgjl(x,t),
b j k ( x , t ) = 2 π ∑ l , s = 1 m K l s g j l ( x , t ) g k s ( x , t ) b_{jk}(\mathbf{x},t)=2\pi\sum_{l,s=1}^{m}K_{ls}g_{jl}(\mathbf{x},t)g_{ks}(\mathbf{x},t) bjk(x,t)=2πl,s=1mKlsgjl(x,t)gks(x,t)

由此可以得出, f 1 = x 2 , f 2 = − h ( x 1 , x 2 ) f_{1}=x_{2},f_{2}=-h(x_{1},x_{2}) f1=x2,f2=h(x1,x2), g 1 j = 0 ( j = 1 , 2 , 3 ) g_{1j}=0(j=1,2,3) g1j=0(j=1,2,3)
g 21 = x 1 , g 22 = x 2 , g 23 = 1 g_{21}=x_{1},g_{22}=x_{2},g_{23}=1 g21=x1,g22=x2,g23=1, n = 2 , m = 3 n=2,m=3 n=2,m=3.
则一、二阶导数矩:
a 1 = x 2 , a 2 = − h ( x 1 , x 2 ) + π K 22 x 2 a_{1}=x_{2},a_{2}=-h(x_{1},x_{2})+\pi K_{22}x_{2} a1=x2,a2=h(x1,x2)+πK22x2
b 11 = 0 , b 12 = 0 , b 21 = 0. b 22 = 2 π K 11 x 1 2 + 2 π K 22 x 2 2 + 2 π K 33 . b_{11}=0,b_{12}=0,b_{21}=0.b_{22}=2\pi K_{11}x_{1}^{2}+2\pi K_{22}x_{2}^{2}+2\pi K_{33}. b11=0,b12=0,b21=0.b22=2πK11x12+2πK22x22+2πK33.
从而得到FPK方程:
∂ ∂ t p + ∂ ∂ x 1 ( x 2 p ) + ∂ ∂ x 2 { [ ( − h ( x 1 , x 2 + π K 22 x 2 ] p } − π ∂ 2 ∂ x 2 2 [ ( K 11 x 1 2 + K 22 x 2 2 + K 33 ) p ] = 0 \frac{\partial}{\partial t}p+\frac{\partial}{\partial x_{1}}(x_{2}p)+\frac{\partial}{\partial x_{2}}\left \{ [(-h(x_{1},x_{2}+\pi K_{22}x_{2}]p \right \}-\pi \frac{\partial ^{2}}{\partial x_{2}^{2}}[(K_{11}x_{1}^{2}+K_{22}x_{2}^{2}+K_{33})p]=0 tp+x1(x2p)+x2{[(h(x1,x2+πK22x2]p}πx222[(K11x12+K22x22+K33)p]=0.

相应的也可以得到Stratonovich或Ito方程,它们得到的FPK方程都和上式一致。

在这里插入图片描述
在这里插入图片描述

3.绘制随机微分方程概率密度曲线

在这里插入图片描述

%(b)
clc;clear;
% 定义需要绘制的变量和参数
t_values = linspace(0.01, 5, 100); % 在t轴上均匀采样100个点,保证不会出现除数为0的情况
x_values = linspace(-5, 5, 200); % 在x轴上均匀采样200个点
[x_mesh, t_mesh] = meshgrid(x_values, t_values); % 创建网格

% 计算概率密度函数
p = (1./sqrt(2*pi.*t_mesh)).*cosh(x_mesh).*exp(-0.5./t_mesh).*exp(-(x_mesh.^2)./(2*t_mesh));

% 绘制演化图
mesh(x_mesh, t_mesh, p); % 使用surf函数绘制三维曲面图
xlabel('x'); ylabel('t'); zlabel('p(x,t)'); % 添加轴标签
title('Probability Density Evolution'); % 添加标题

在这里插入图片描述

在这里插入图片描述
未完待续…

本文来自互联网用户投稿,该文观点仅代表作者本人,不代表本站立场。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如若转载,请注明出处:http://www.coloradmin.cn/o/511969.html

如若内容造成侵权/违法违规/事实不符,请联系多彩编程网进行投诉反馈,一经查实,立即删除!

相关文章

[医学分割比赛] ISBI2023 APIS多模态医学分割比赛总结 + top3解决方案

ISBI2023 APIS多模态医学分割比赛总结 top3解决方案 0.比赛背景1.比赛任务及结果2.第三名方案 - 龙盈智达(北京)科技有限公司(0) Data Preprocessing(1) Data Augmentation(2) Approach(Model)(3) Approach(Data Sampling)(4) Ap…

【QT】学习课-pushButton的使用(1)!

Qt 是一个1991年由Qt Company开发的跨平台C图形用户界面应用程序开发框架。它既可以开发GUI程序,也可用于开发非GUI程序,比如控制台工具和服务器。Qt是面向对象的框架,使用特殊的代码生成扩展(称为元对象编译器(Meta Object Compi…

《WebGIS快速开发教程》写好啦

告诉大家一个好消息,经过我没日没夜,呕心沥血的创作,这本叫做《WebGIS快速开发教程》的书籍终于写好了。这本书适用于还未毕业的学生、以及正在从事传统前后端开发但是想转到WebGIS开发的人。 这本书的特点突出一个“快”和“轻”&#xff0c…

三子棋小游戏---(C语言)

目录 前言: 1.菜单的打印 2.三子棋构思 3.实现三子棋 3.1使用宏的方式定义数组 3.2打印棋盘 3.3玩家下棋 3.4电脑随机下棋 3.5判断结局 ❤博主CSDN:啊苏要学习 ▶专栏分类:C语言◀ C语言的学习,是为我们今后学习其它语言打好基础&am…

Kyligence Zen产品体验——一站式指标平台泰酷辣~

文章目录 一、前言二、为什么需要指标化平台三、什么是Kyligence Zen四、Kyligence Zen新特性五、Kyligence Zen注册篇六、Kyligence Zen体验篇七、Kyligence Zen实战篇7.1 导入数据7.2 创建指标7.3 指标分析 八、Kyligence Zen总结篇九、参考资料 一、前言 随着互联网和物联网…

tomcat集群下的session共享和负载均衡(redis实现)

环境 操作系统:windows tomcat1:Apache Tomcat/7.0.52(8085) tomcat2:Apache Tomcat/7.0.52(8086) jre:1.7.0_80 nginx:nginx-1.20.1(8070) redis…

基于 SpringBoot+WebSocket 无DB实现在线聊天室(附源码)

文章目录 基于 SpringBootWebSocket 无DB实现在线聊天室0 项目说明0.1 样例展示0.2 源码地址 1 WebSocket 简介1.1 HTTP1.2 WebSocket1.2.1 WebSocket 协议1.2.2 WebSocket 交互 2 使用教程2.1 客户端(浏览器)2.1.1 WebSocket 对象2.1.2 WebSocket 事件2…

重装系统后,qt5.11.3升级到qt5.12.6所遇到的问题

前提:重装了系统: c/qt windows10 语音模块TTS异常,数据库缺少驱动 一:语音模块不能播放 qt使用语音模块时,在初始化时出现异常: onecore\com\combase\dcomrem\resolver.cxx(2299)\combase.dll!00007FF8…

Oracle存储过程~封神之路

简介 Oracle 存储过程是 Oracle 数据库中的一种数据处理对象,它可以在数据库中定义一组预定义的 SQL 语句,用于完成特定的数据库操作。存储过程可以被授权的用户调用,并且可以执行多个语句,这些语句可以被视为一个单独的操作&…

“深圳首届十大金口碑人物”优必选科技创始人兼CEO周剑获此殊荣

深圳晚报社联合深圳市诚商信用评级有限公司、深圳市诚信营商促进会和中国善网,共同举办了首届“金口碑”评选活动。活动涵盖多个领域,历经多个环节的评定和实地走访,最终有10名个人、20家企业和70家商户成功获得“深圳首届十大金口碑人物”、…

Visual C++实现推箱子游戏的核心算法设计与实现(附源码和和资源)

需要源码和资源请点赞关注收藏后评论区留言私信~~~ 在前面的博客中已经讲解了推箱子游戏的菜单和各种对话框的实现,下面对推箱子游戏的核心算法设计和实现进行讲解 一、地图文件读取模块的设计与实现 地图文件读取模块,主要负责将地图文件进行读取&…

【AI大模型】“讯飞星火”认知大模型正式发布 计划10月底赶超ChatGPT

文章目录 前言你使用过这种对话式AI吗?有什么看法或感受?“讯飞星火大模型将超越chatgpt?”这类型的人工智能对现在的社会有什么意义?这类型的人工智能,未来前景如何?申请体验写在最后 前言 5月6日&#xf…

科普:跨链桥是如何被黑的?

科普:跨链桥是如何被黑的? [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-hbvPVTkb-1666679410875)(htTPS://tva1.sinaimg.cn/large/e6c9d24ely1h4r0auxvmxg20tr04ojug.gif)] 跨链桥是一种允许两个独立区块链网络之间进行通信…

【git】用好 stash,工作超nice

一、介绍 如果修改后的内容还不想commit,就可以用git stash命令。它会将工作区和暂存区中的修改(也就是还没commit的内容)都会被保存到堆栈里,并在之后恢复到任意指定的分支上。 二、应用场景 1、在分支a进行开发feature 1时,突然需要紧急…

【AI大模型】讯飞星火大模型能否超越chatgpt?

文章目录 前言你使用过这种对话式AI吗?有什么看法或感受?“讯飞星火大模型将超越chatgpt?”这类型的人工智能对现在的社会有什么意义?这类型的人工智能,未来前景如何?申请体验写在最后 前言 5月6日&#xf…

机器学习-10 聚类算法

聚类算法 算法概括聚类(clustering)聚类的概念聚类的要求聚类与分类的区别 常见算法分类聚类算法中存在的问题 距离度量闵可夫斯基距离欧式距离(欧几里得距离)曼哈顿距离切比雪夫距离皮尔逊相关系数余弦相似度杰卡德相似系数 划分…

建造者模式详解:建造随意搭配的肯德基套餐

一、简介 建造者模式(Builder Pattern)是五种创建型设计模式之一,它将一个复杂对象的构建与其表示分离,使得同样的构建过程可以创建不同的表示。这句话怎么理解呢:一个对象的构建过程本质上就是这个对象包含的所有成员…

一觉醒来IDEA感觉不香了,AI智能编程工具Cursor使用

一、简介 为使用人工智能编程而构建的编辑器,一款人工智能编程软件、智能Ai代码生成工具。 它有什么特点呢? 集成了GPT-4,国内可用,不仅有ChatGPT的聊天功能,还有强大的自动代码生成能力,简直是编码神器。 …

Java中常见的几种数组排序方法

这篇文章总结一下我学习到的几种常见的数组排序方法 冒泡排序 冒泡排序在我看来是最简单、最基本的排序方法,我们应当将冒泡排序的原理和代码熟记于心。 冒泡排序的原理十分简单:用数组的第一个元素和第二个元素进行比较,将大的放到后面&a…

【Java编程系列】Minio实现文件上传下载

热门系列: 【Java编程系列】Amazon S3实现文件上传下载 目录 热门系列: 1、前言 2、Minio实战代码 2.1 Minio环境部署 2.2 Minio的Sdk对接实现 2.2.1 Minio Maven依赖 2.2.2 minio配置与初始化 2.2.3 上传文件 2.2.4 下载文件 2.2.5 生成文件…