一个常规的微分方程。

$\frac { d y } { d t } = a y + e ^ { s t }$ $y = y(0) at t = 0$ Look for $y_p = Ye^{st}$

把$\frac { d y } { d t }$比做银行的钱的增长率，$y$就是银行的利息，$e ^ { s t }$代表挣到的钱，是指数形式的，$s$代表souce，源，t代表时间。

带入$y_p = Ye^{st}$，得到$Y = \frac {1}{s-a}$

则 $y(t) = \frac{e^{st}}{s-a} + null solution$

null solution的意思是，输入项没有了，也就是$e^{st} = 0$，求$\frac { d y } { d t } = a y$，得$y = e^{at}$

所以，通解

$y(t) = \frac{e^{st}}{s-a} + Ce^{at} = y_p + y_n\tag{1}$

y_p意思是y_particular，y_n意思是y_null。

将$y = y(0) at t = 0$带入得特解

$y(0) = \frac {1}{s-a} + C$

则$C = y(0) - \frac {1}{s-a}$，带回（1）得，将y(0)单独提出来

$\begin{aligned} y ( t ) &\left. = \frac { e ^ { s t } } { s - a } + \left[ y ( 0 \right) - \frac { 1 } { s - a } \right] e ^ { a t } \ & = \frac { e ^ { s t - e ^ { a t } } } { s a } + y ( 0 ) e ^ { a t } \ & = y_{vp} + y_n \end{aligned}$

y_vp是y_very particular。

特接并不唯一。

这个null solution很好，可以看做是从初始状态开始的一个变化，并包含初始状态信息。

这个特别的解很好，有这样的性质，当t=0时，特解同时=0，这样就得到了通解。

但是有一个例外，当s=a时，分母为0，同时分子也为0，这时，就出现了震荡（Resonance）。这时候需要用到洛必达（L‘ Hopital）法则。

IF s = a，then

$y_{vp}+y_n = t e^{at} + y(0)e^{at}\tag{2}$

由此得到公式(2)中的特解部分。

test

如何给人分享一个东西让人家觉得珍贵，而自己又付出的较少。

利用迭代效应影响他人。

比如茶叶，普洱茶流传着一句至理名言，“越陈越香，越老越贵”。说的是茶存放的时间越长它的品质越好，相应的茶会升值。

利用这个效应来送礼就很合适，送的可以有时间的加成。

朋友请自己喝百茶堂的老茶，我们心里觉得这茶就很贵，很好，自然入喉后仔细的品味一番，然后就会不自然的说这是好茶。

而实际上，朋友买的时候的价格并没有特别贵重，这就是时间的味道。

反过来想问题就是，现在就要为未来请客积攒礼物，如何挑选越陈越值钱的礼物是很关键的，这个需要再仔细思量。

喝到一款很好喝的茶，而并不是百茶堂，忽有感，百茶堂的茶也并没有那么优秀，而在王晨那里喝到很好的茶，年份也好，但就不会觉得特别优秀，这就是如何营造时间的意义方面的话题了。

关于克服拖延症的5个建议

1. 改变想法

   一个事情我今天不去做，就不会再去做了，想到一个事情，若必须去做，那么今天就去做，可以做任务分解，然后分解到每一天，若任务太多，将任务分级，做更重要的。

2. 每天进行任务划分

   每天做任务计划，然后按照计划去完成，每天都要完成，最开始可以划的粗糙一点，任务量轻一点，然后逐步加大。

3. 仪式感

   每天做学习的的任务之前要给自己定一个仪式，就和睡眠仪式一样的，有助于进入状态

4. 每日打卡

   对每天的完成情况打卡，有助于积累成就感。

5. 对自己加油！

   加油吧！

The Calculus You Need

学习微分方程需要用到一些前置的知识，主要有以下几方面

1. 特殊的一些函数的微分

$x^n$ $sin(x)$ $cos(x)$ $e^x$ $ln(x)$

清楚以上特殊函数的微分。

2. 微分操作定律

$f+g$ $f(x)+g(x)$ $\frac { f ( x ) } { g ( x ) }$ $f(g(x))$

清楚以上方程式的微分。

3. 一个求解公式

用$y ( t ) = \int _ { 0 } ^ { t } e ^ { t - s } 2 ( s ) d s$用来求解$\frac { d } { d x } \int _ { 0 } ^ { x } y ( t ) d t = y ( x )$

我们可以对其证明，对左边的式子左右两端进行微分。

$\frac { d y ( t ) } { d t } = \frac { d } { d t } \left( e ^ { t } \cdot \int _ { 0 } ^ { t } e ^ { - s } q ( u ) d s \right)$

得到

$\frac { d y ( t ) } { d t } = e ^ { t } \int _ { 0 } ^ { t } e ^ { - s } g ( x ) d s + g(t)$

即

$\frac { d } { d x } \int _ { 0 } ^ { x } y ( t ) d t = y ( x )$

4. Taylor级数

我们用一阶导数来形容斜率，我们用二阶导数来形容bending，我们用三阶导数以及更高阶导数去修正曲线，在$t+\Delta t$这段，也就是$\frac{df}{dt}$段更加拟合原函数$f(x+\Delta t)$

$f(t+\Delta t) = f(t) + \Delta t \frac {df(t)}{dt} + \frac{1}{2}(\Delta t)^2 \frac{d^{2}t}{dt^{2}}+ \cdots + \frac { 1 } { n ! } ( \Delta t ) ^ { n } \frac { d ^ { n } f } { d t ^ { n } } + \cdots$

Overview of Differential Equations

https://www.youtube.com/watch?v=ghjOS7Q82s0&list=PLUl4u3cNGP63oTpyxCMLKt_JmB0WtSZfG&index=2

一个典型的一阶的微分方程

$\frac { d y } { d t } = a y + q ( t )$

该方程描述的是y的变化率由y自己的解来决定，同时会在等式右端引入一个激励。这是一个线性的y的方程，因为y是一次的。

$\frac { d y } { d t } = f ( y )$

这可能是一个非线性的，因为$f(y)$可以是y的任何函数，可以是$y^2$、$sin(y)$、$e^y$等等可能形式。

这里可以形象的理解，y是弹簧的长度，$\frac { d y } { d t }$是弹簧长度的变化率。

二阶微分方程

$\frac { d ^ { 2 } y } { d t ^ { 2 } } = - k y$

二阶微分方程可理解为加速度，也是区间的弯折率（beding of a curve），一阶微分方程给的是曲线的方向，是上升还是下降，二阶微分方程给的是曲线向上弯还是向下弯，也就是凸函数还是凹函数。也可以理解为力的变化加速度，而线性的与y有关。

更一般的形式

$m y ^ { \prime \prime } + b y ^ { \prime } + k y = f(t)$

$m y ^ { \prime \prime }$理解为加速度乘以质量，第二项$b y ^ { \prime }$表示阻尼，速度的变化量，第三项是一些强制条款，取决于y自己，$f(t)$是一些外部因素。

我们可解的是线性常微分方程，Linear Constant Coefficients，这样的微分方程我们认为是good equations。

对于不好的方程，比如非线性方程，比如不是常微分方程，即系数（m/k/b）是不变的方程，则需要用数值方法求解。

一个系统经常不止是一个方程，是Systems of n equations.

$\frac { d y } { d t }= A y$其中A是矩阵 $\frac { d ^ { 2 } y } { d t ^ { 2 } } = - S y$

特征值和特征向量。

数值解法很重要Numerical Solutions

目前很好用的方法是Ode45.

而Euler方法可看做是Ode1，Ode45具有更高的精度。“Ode 4 and 5”

偏微分方程，有2个变量，以下是热力学方程，时间以及x是空间向量。

$\frac { \partial u } { \partial t } = \frac { \partial ^ { 2 } u } { \partial x ^ { 2 } }$

以上是一个很重要常微分方程。（PDE，Partial differential equations）偏微分方程

$\frac { \partial ^ { 2 } u } { \partial t ^ { 2 } } = \frac { \partial ^ { 2 } u } { \partial x ^ { 2 } }$

上面是波动方程。

拉普拉斯方程。

$\frac { \partial ^ { 2 } u } { d x ^ { 2 } } + \frac { \partial ^ { 2 } u } { d y ^ { 2 } } = 0$

对于Cassie电弧方程的离散化。

引用

自己编写了一个s函数，有几个参数引用了自己写的几个函数，在脚本中可以正确运行，但在写成s函数，进行 simulink 仿真的时候，已知提示“too many input auguments”,不知道怎么回事。

经过调试发现，把那几个参数换成常值，simulink 可以正常运行。说明：s函数中，不能直接调用自己写的函数。

那么问题来了：

在用 sfuntmpl 模板用m语言写的s函数中，不能直接调用自己写的函数，怎么办？

具体的解决办法是：

1、打开主页上的‘set path’选项，把自己写的函数所在的文件夹加到搜素路径中。

2、在matlab的安装文件路径 X:\Program Files\MATLAB\R2012\toolbox\ 下，新建一个文件夹 ‘myTool’。把自己写的函数复制到此文件夹下，然后把次文件夹加到路径中。

任选一种方法即可。个人比较推荐2。

在进行电弧模型建模时候，传统的cassie电弧模型的写法都是连续的，但是连续状态的采样率设置上就出现了问题，无法按照自己期望的采样率进行仿真，所以需要对模型进行离散化处理。

转自csdn

建立微分方程只是解决问题的第一步通常需要求出方程的解来说明实际现象，并加以检验。如果能得到即系形式的解固然是便于分析和应用的，但是我们知道，只有线性常系数微分方程，并且自由项是某些特殊类型的函数时，才可以肯定得到这样的解，而绝大多数变系数方程、非线性方程都是所谓“解不出来”的，即使看起来非常简单的方程如$\frac{d y}{d x}=y^{2}+x^{2}$，对于用微分方程解决实际问题来说，数值解法就是一个十分重要的手段。

1. 常微分方程离散化

下面主要讨论一阶常微分方程的初值问题，其一般形式是

在下面讨论中，我们总假定函数$f(x,y)$连续，且关于y满足Lipschitz条件，即存在常数L，使得

这样，由常微分方程理论知，处置问题（1）的解必定存在唯一。

数值解法

所谓数值解法，就是求问题（1）的解$y(x)$在若干点

$a=x_{0}<x_{1}<x_{2}<\cdots<x_{N}=b$处的近似值$y_n(n=1,2,…,N)$，$y_n(n=1,2,…,N)$称为问题（1）的数值解，$h_n=x_{n+1}-x_n$称为由$x_n$到$x_{n+1}$的步长，今后如无特别说明，我们总取步长为常量h。

建立数值解法，首先要将微分方程离散化，一般采用以下几种方法：

（i）用差商近似导数——缠粉方程初值问题

若用向前差商$\frac{y\left(x_{n+1}\right)-y\left(x_{n}\right)}{h}$代替$y^{\prime}\left(x_{n}\right)$带入（1）中的微分方程，则得

简化得

如果用$y(x_n)$的近似值$y_n$带入上式右端，所得结果作为$y(x_{n+1})$的近似值，记为$y_{n+1}$，则有

这样，问题（1）的近似解可通过求解下述问题

按式（3）由初值$y_0$可逐次算出$y_1,y_2…$。式（3）是个离散化的问题，成为差分方程初值问题。

（ii）用数值积分方法

将问题（1）的解表成积分形式，用数值积分方法离散化。例如，对微分方程两端积分，得

右边的积分用矩形公式或梯形公式计算。

（iii）Taylor多项式近似

将函数$y(x)$在$x_n$处展开，取一次Taylor多项式近似，则得

再将$y(x_n)$的近似值$y_n$代入上式右端，所得结果作为$y(x_{n+1})$的近似值$y_{n+1}$，得到离散化的计算公式

以上三种方法都是将微分方程离散化的常用方法，每一类方法又可导出不同形式的计算公式。其中Taylor展开法，不仅可以得到求数值解的公式，而且容易估计截断误差。

2 欧拉（Euler）方法

2.1 向前Euler公式、向后Euler公式

Euler方法就是用查分方程初值问题（3）的解来近似微分方程初值问题（1）的解，即由公式（3）依次算出$y(x_n)$的近似值$y_n$，（$n=1,2,…$）。这组公式求问题（1）的数值解成为向前Euler公式。

如果在微分方程离散化时，用向后差商代替导数，即$y^{\prime}\left(x_{n+1}\right) \approx \frac{y\left(x_{n+1}\right)-y\left(x_{n}\right)}{h}$，则得计算公式

用这组公式求问题（1）的数值解成为向后Euler公式。

向后Euler法与Euler法形式上相似，但实际计算时却复杂的多。向前Euler公式是显式的，可直接求解。向后Euler公式的右端含有$y_{n+1}$，因此是隐式公式，一般要用迭代法求解，迭代公式通常为

2.2 Euler方法的误差估计

对于向前Euler公式（3）我们看到，当n = 1,2,……时公式右端的$y_n$都是近似的，所以用它计算的$y_{n+1}$会有累积误差，分析累积误差比较复杂，这里先讨论比较简单的所谓局部截断误差。

假定用（3）式右端的$y_n$没有误差，即$y_n=y(x_n)$，那么由此算出

局部截断误差指的是，按（7）式计算由$x_n$到$x_{n+1}$这一步的计算值$y_{n+1}$与精确值$y(x_{n+1})$之差$y(x_{n+1})-y_{n+1}$。为了估计它，由Taylor展开得到的精确值$y(x_{n+1})$是

(7)、(8)两式相减（注意到$\left.y^{\prime}=f(x, y)\right)$）得

即局部截断误差是$h^2$阶的，而数值算法的精度定义为：

若一种算法的局部截断误差为$O(h^{p+1})$，则称该算法具有p阶精度。

显然p越大，方法的精度越高。式（9）说明，向前Euler方法是一阶方法，因此它的精度不高。

3 改进的Euler方法

3.1 梯形公式

利用数值积分方法将微分方程离散化时，若用梯形公式计算式（4）中之右端积分，即。

这就是求解初值问题（1）的梯形公式。

直观上容易看出，用梯形公式计算数值积分要比矩形公式好。梯形公式为二阶方法。 梯形公式也是隐式格式，一般需用迭代法求解，迭代公式为

如果实际计算时精度要求不太高，用公式（10）求解时，每步可以只迭代一次，由此导 出一种新的方法—改进 Euler 法。

3.2 改进Euler法

按式（5）计算问题（1）的数值解，如果每步只迭代一次，相当于将Euler公式与梯形公式结合使用：先用Euler公式求$y_n+1$的一个初步近似值$\bar{y}_{n+1}$，称为预测值，然后用梯形公式校正求得近似值$y_n+1$，即

式（11）称为由Euler公式和梯形公式得到的预测——校正系统，也叫改进Euler法。

为了便于编制程序上机，式（11）常改写成

改进Euler法是二阶方法。

4 龙格——库塔（Runge-Kutta）方法

这部分以后再补，暂时用不上。

https://blog.csdn.net/qq_29831163/article/details/89703598

5 Matlab解法

https://blog.csdn.net/qq_29831163/article/details/89703911

5.1 非刚性常微分方程的解法

Matlab 的工具箱提供了几个解非刚性常微分方程的功能函数，如 ode45，ode23， ode113，其中 ode45 采用四五阶 RK 方法，是解非刚性常微分方程的首选方法，ode23 采用二三阶 RK 方法，ode113 采用的是多步法，效率一般比 ode45 高。 Matlab 的工具箱中没有 Euler 方法的功能函数。

（1）对简单的一阶方程的初值问题

自己编写改进的Euler方法函数eulerpro.m如下：

function [x,y]=eulerpro(fun,x0,xfinal,y0,n);
if nargin<5,n=50;end 
h=(xfinal-x0)/n;
x(1)=x0;y(1)=y0;
for i=1:n
    x(i+1)=x(i)+h;
    y1=y(i)+h*feval(fun,x(i),y(i));
    y2=y(i)+h*feval(fun,x(i+1),y1);
    y(i+1)=(y1+y2)/2;
end 

1. 可能提示没有对应的function

看一看对应的function在setup模块中是否被引用。

1. 可能提示Dwork不能有空的。。。

看一下Dwork设定参数的个数，与实际设定的Dwork参数个数是否相符。

1. 不同的function中的参数无法互相通用，也是程序无法正确使用的原因。

Dwork有效性试错

如果Dwork可以看做是全局变量，那么在InitializeConditions中的Dwork的值是可以带入到Output子函数的，现在试试。

设置0.1s，step应该变值，但是很显然没有变化，所以这样是不行的。

也就是说InitializeConditions并不是执行一次就初始化一次么？

调整将Dwork放到Derivatives模块中，这样就相当于原来的了，应该就好用了。

function Derivatives(block)

tau = block.DialogPrm(1).Data;
Uc = block.DialogPrm(2).Data; 

if block.CurrentTime <= block.DialogPrm(3).Data
    block.Dwork(3).Data = 0;
else
    block.Dwork(3).Data = 1;
end
    
block.Derivatives.Data = block.Dwork(3).Data * 1/tau * (block.InputPort(1).data^2 / Uc^2 - 1);
%endfunction

这样是可以的。