Hamilton–Jacobi 方程特征线法求精确解
与 Burgers 方程为代表的守恒律方程类似,Hamilton–Jacobi 方程也可以利用特征线方法进行精确求解。 一维 Hamilton–Jacobi 方程如下 \[ \left\{ \begin{aligned} u_t + H\left(u_x\right) &= 0\\ u(x,0) &= u_0(x) \end{aligned} \right. \] 其中定义域 \(x \in \mathbb{R}\),通常假定初值有周期性,精确解此时也具有周期性。为了简化记号,引入符号 \(p = u_x\)。 我们将考虑两种具体情况: \(H(p) = \frac12 (p+\alpha)^2\),其中 \(\alpha > 0\) 为常数,此时 \(H(p)\) 是一个严格凸函数; \(H(p) = - \cos(p + \alpha)\),此时 \(H(p)\) 是一个非凸非凹的函数。 特征线方程 从\((\xi,0)\)发出的曲线记作 \(l_\xi\),曲线参数为 \(s\),方程为 \(x = x(s)\),其中...
Numpy 学习笔记之三
NumPy 提供了比较丰富的线性代数运算功能,主要集中在 numpy.linalg 模块中,下面罗列一些基础的操作,主要针对一维和二维数组,对于其它高维数组的语义可能不太符合直觉。 重要函数 点积 np.dot np.dot() 函数可以计算两个向量的点积,最基本的用法例如 1234a = np.array([1, 2, 3], dtype=np.float64)b = np.array([4, 5, 6], dtype=np.float64)np.dot(a, b) # np.float64(32.0)# 1*4 + 2*5 + 3*6 = 32 np.dot(a, b)的完整语义需要根据两个数组的维度来区分: 第一类情况: 如果其中一个是标量,那么就是进行(逐元素的)乘法,等效于 np.multiply(a, b) 或 a * b。 第二类情况: 如果 a 是 \(N\) 维数组,b 是一维数组,那么会沿着两者的最后一个轴求和,得到 \(N-1\) 维数组。 (特例,标准情况)如果两者都是尺寸为...
Numpy 学习笔记之二
数组的拷贝 Python 本身有深拷贝和浅拷贝等复杂的概念,对于多层列表等数据结构,对应的行为差异非常明显。 但是 Numpy 的 ndarray 数组无论形状如何,在内存中始终是连续存储的,高维数组并不是数组的数组, 因此并没有深拷贝和浅拷贝的明显差异。 和其它 Python 数据一样,在函数传参过程中对 ndarray 数组的传递并不会触发任何的拷贝行为。 为了获取一个 ndarray 数组的完整拷贝,可以使用 ndarray.copy() 方法,例如 12345a = np.array([1, 2, 3])b = a.copy()b[0] = 100a # array([1, 2, 3]) 上述测试代码表明,拷贝得到的是完全独立的数组,修改不会相互影响。 也可以使用numpy.copy()函数,例如 12a = np.array([1, 2, 3])b = np.copy(a) ndarray.copy() 方法和...
Numpy 学习笔记之一
Numpy 是一个用于科学计算的 Python 扩展库,虽然不是 Python 标准库,但是已经是 Python 在科学计算、机器学习等领域毋庸置疑的基础库,很多主流的 Python 库都依赖 Numpy。 1import numpy as np 一个循序渐进的入门笔记的内容组织比较麻烦,我也并不打算这么做,这个笔记的内容会围绕一些关键点展开, 主要参考 Numpy 的官方资料。 ndarray 基本概念 np.ndarray (别名 np.array)是 Numpy 的核心数据类型,可以用于存储 n 维数组,并提供了大量相关操作方法。 它与 Python 的 list 类型类似,但是存在更多限制:只能存储同类型数据,并且尺寸是固定且规则的,即不允许各个元素长短不一的数组。 这些限制使得 ndarray 的运算非常高效(直接调用底层的 C/C++ 实现)。 但是有一点不足的是:Numpy 没有提供稀疏矩阵的数据结构,需要时可以使用其它包提供的稀疏矩阵。 我们主要关注一维数组和二维数组,例如 123array([1, 2, 3]) # shape=(3,)array([[1,...
MATLAB 高性能编程笔记
记录一些常用的MATLAB编程技巧/规范,目标是写出高效并且可维护的MATLAB代码。主要参考官方文档中的提升性能的方法。 这里的讨论只关注运算效率,假设内存总是足够的,并且不涉及并行计算和GPU。 关于代码结构 我们有很多方式执行代码:命令行 vs 脚本 vs 函数,虽然绝大多数代码在不同方式中执行都是等效的,但是考虑优化就不是一回事了,通常的运行效率关系为:命令行 < 脚本 < 函数,也就是说函数的执行速度通常是最快的。 基于模块化编程的思维,避免使用过大的单一文件,考虑将其中重复使用的功能拆分为简洁的函数文件,这种做法可以降低首次运行的成本。 优先使用局部函数而非嵌套函数,嵌套函数可以直接访问外层函数的变量(按照引用捕获),这会影响效率,更好的做法是将所有需要的变量以函数参数的形式显式传递。 关于数组操作 预分配+向量化 两个核心原则:预分配 + 向量化 预分配数组:提前分配数组大小,避免在循环中动态扩展数组。 向量化计算:尽可能使用矩阵和向量运算,避免 for 循环。 我们可以对 for 循环中的几种常见写法进行简单的对比:...
Python 理解切片的底层原理
在学习 Python 的字符串和列表所支持的切片操作,以及 Numpy 所支持的各种花式索引操作时,总是会感到非常困惑, 究其原因,就是没有理解这背后的底层原理,因此有必要专门学习一下。 __getitem__ 和 __setitem__ 方法 Python 为类型的 [] 操作提供支持的方法是 __getitem__ 方法和 __setitem__ 方法,分别提供读写的功能。 例如我们可以实现一个支持 [] 读操作的实验类型: 1234class Demo: def __getitem__(self, index): print(f"type(index) = {type(index)}") print(f"index = {index}") 尝试一下基本的索引操作 12345678910111213a = Demo()a[0]# type(index) = <class 'int'># index = 0a[-1]#...
MATLAB 单元测试
学习一下关于 MATLAB 单元测试的内容,为了维护一份健壮的代码库,单元测试是必不可少的,参考官方文档。 MATLAB 提供的测试主要包括三种风格: 基于脚本 基于函数 基于类 基于脚本的单元测试 可以通过一个脚本对指定功能进行测试,脚本名称必须以 test 开头或结尾,不区分大小写,否则测试文件可能被忽略。 每一代码节(%%)作为一个测试单元,随后的文本视作测试单元的名称,否则MATLAB会提供一个默认名称。 如果脚本中不含代码节,那么整个文件会被视作一个测试单元,测试单元的名称为脚本名称。 在测试单元中主要通过 assert...
MATLAB 帮助系统
这篇笔记只为了解决两个问题:如何基于命令行使用MATLAB的帮助文档系统?如何让自己编写的代码(函数和类)适配MATLAB的帮助系统? 获取帮助 doc & docsearch 打开 MATLAB 文档浏览器(最新版已经改成使用系统浏览器),显示帮助中心对应的完整页面。 打开帮助中心首页 1doc 打开abs对应的页面 1doc abs 在帮助中心搜索关键字 1docsearch abs 这两个命令都是完全基于离线版的MATLAB帮助中心,和在线的MATLAB帮助中心几乎一致(除了版本和语言可能不同),通常不支持对用户自定义代码的帮助。 help help 命令可以在命令行中显示简要帮助信息。 help name 会显示 name 指定的功能的帮助文本,例如函数、方法、类、工具箱、变量或命名空间。 help 则会显示与先前操作相关的内容。 例如查找函数的帮助信息 1help abs 输出内容如下 123456789abs Absolute value. abs(X) is the absolute value of the elements of X....
MATLAB 踩坑记录
记录一下 MATLAB 踩过的小坑,MATLAB的各种内置函数的用法实在是太奇怪了,一个小细节就可以搞出各种问题。 矩阵元素个数 除了最常见的size(A)可以获取数组的尺寸,还有如下函数: numel(A) 可以获取(任意维度数组)所有元素个数; length(A) 可以获取向量的长度(也就是元素个数),但是扩展到高维数组的行为是非常反直觉的——返回最大维数长度,相当于max(size(A))。 矩阵偏移 函数 circshift 可以用于向量或矩阵的循环偏移,不会修改原始数据。 向量偏移 123456A = [1, 2, 3, 4, 5];circshift(A, 1)% [5, 1, 2, 3, 4]circshift(A, -1)% [2, 3, 4, 5, 1] 矩阵的偏移就比较复杂了,按照前面的做法会以行为整体进行偏移!!! 1234567A = [1 2 3; 4 5 6; 7 8 9];circshift(A, 1);% 7 8 9% 1 2 3% 4 5 ...
Cpp 多态方案对比——回字的四种写法
虽然C++没有直接提供interface,但是却提供了虚函数、模板类型等各种语法,使得我们可以用各种方式实现多态,这里我们不区分动态多态和静态多态,而是从设计一个框架的角度,分别使用四种方案实现: 虚函数(最简单直接的方式) std::function CRTP(最晦涩的方式) deducing this(可以视作CRTP的简化,不再需要将派生类作为模板参数传递,要求C++23) 需求 我们考虑这样一个需求: 基类A包括:(不可实例化) 主方法run:调用func1,func2和func3 实现方法func1 实现方法func2(多态,允许子类修改) 实现方法func3(多态,子类必须实现) 派生类B1:继承A 实现方法func3(多态,允许子类修改) 派生类B2:继承A 实现方法func3(多态,允许子类修改):调用func4 实现方法func4(多态,允许子类修改) 具体类C1:继承B2 实现方法func4 具体类C2:继承B2 实现方法func3 最终我们直接通过各种对象自身来调用run方法,达到如下效果: 12345678Running B1...