精华内容
下载资源
问答
  • 针对低矮小开口剪切型剪力墙,研究...最后应用Ritz变换,将低矮小开口剪力墙的自振特性分析问题转化为一组非线性代数方程的求解问题,进而求得该类剪力墙自振特性的近似半解析解。算例结果表明,该方法对以剪切变形影响为
  • 电渗流在粗糙板上的流动特性分析。通过微扰法求得粗糙板电渗流半解析解,分析结果。
  • 基于已有的地基反力系数为常数时桩的变形和内力的解析解,以及桩的响应的连续性,提出了半解析解。采用Fortran语言编制了计算程序,并通过特例情况下的解与已有解进行对比来验证所得解的可靠性。计算结果表明:桩长...
  • 首先对Fredlund的非饱和土一维固结理论进行简化,由得到的...采用Crump方法编制程序实现Laplace 逆转换,得到了时间域内的超孔隙水压力、超孔隙气压力、土层沉降的半解析解;引用典型算例,对单面排水排气情况,与已有的解析
  • 从三维轴对称土模型出发,考虑桩周土体的三维波动效应,对均质滞回材料阻尼土中完整端承桩在垂直协和激振力作用下的纵向振动特性进行分析,求解得到桩顶...求得半正弦脉冲激振力作用下桩顶速度导纳时域响应半解析解
  • 在指数形式渗流退化为达西渗流的条件下,将半解析法的计算结果与解析解进行比较,验证了半解析法计算结果的可靠性。最后,结合某双层地基固结实例对不同参数时的固结性状进行了分析。结果表明:用半解析法计算基于指数...
  • 考虑桩周土体的三维波动效应,分析成层黏性材料阻尼土中桩顶受任意轴向激振力作用时变阻抗桩的纵向振动特性,求得桩顶频域响应解析解、复刚度和速度导纳,分析了底部土层模量变化和桩身阻抗变化对桩纵向振动特性的...
  • 通过采用数值离散与成层地基线性固结解析解相结合的半解析方法得到一系列的固结曲线.分析仅在某一参数发生变化时的一维固结性状.结果表明:固结速率在固结初期随渗流指数的增加而增加,在固结后期随渗流指数的增加...
  • BAW美式期权定价半解析解 参考文献:G. Barone-Adesi & R. E. Whaley, Efficient Analytic Approximation of American Option Values, The Journal of Finance, No 2 (June 1987), 301-320 基本思想:美式期权...

    BAW美式期权定价半解析解

    参考文献:G. Barone-Adesi & R. E. Whaley, Efficient Analytic Approximation of American Option Values, The Journal of Finance, No 2 (June 1987), 301-320

    基本思想:美式期权不同于欧式期权之处在于,美式期权可以在到期日之前随时执行,因而美式期权的价值高于欧式期权。美式期权价值高于欧式期权的部分叫做美式期权溢价 ϵ c \epsilon_c ϵc,如下式所示:
    ϵ c ( S , T ) = C ( S , T ) − c ( S , T ) , \epsilon_c\left( S , T \right) = C \left( S , T \right) -c \left( S , T \right) , ϵc(S,T)=C(S,T)c(S,T),
    其中 C ( S , T ) C \left( S , T \right) C(S,T)表示美式看涨期权的价值; c ( S , T ) c \left( S , T \right) c(S,T)表示欧式看涨期权的价值。

    根据Feynman_Kac定理,在无套利假设条件下,欧式期权价值和美式期权价值满足椭圆的偏微分方程,因而,风险溢价也应该满足以下的相同形式的椭圆偏微分方程:
    1 / 2 σ 2 S 2 ε S S − r ε + b S ε S + ε t = 0 (1) 1 / 2 \sigma^{2} S^{2} \varepsilon_{S S}-r \varepsilon+b S \varepsilon_{S}+\varepsilon_{t}=0\tag{1} 1/2σ2S2εSSrε+bSεS+εt=0(1)
    做如下的变量替换:

    τ = T − t \tau = T -t τ=Tt 代表距离到期日T的时间;

    M = 2 r / σ 2  和  N = 2 b / σ 2 M=2 r / \sigma^{2} \text { 和 } N=2 b / \sigma^{2} M=2r/σ2  N=2b/σ2

    则式(1)可以重写为
    S 2 ε S S − M e + N S e S − ( M / r ) ε τ = 0 (2) S^{2} \varepsilon_{S S}-M_{e}+N S e_{S}-(M / r)\varepsilon_{\tau}=0\tag{2} S2εSSMe+NSeS(M/r)ετ=0(2)
    定义美式期权的提前行权溢价满足如下的形式:
    ϵ c ( S , K ) = K ( τ ) f ( S , K ) \epsilon_c(S, K)=K(\tau) f(S, K) ϵc(S,K)=K(τ)f(S,K)
    其中 K ( T ) = 1 − e − r τ , K \left( T \right) = 1 - e ^ { - r \tau} , K(T)=1erτ, 则(2)式可以重写为
    S 2 f s s + N S f s − ( M / K ) f − ( 1 − K ) M f K = 0 (3) S^{2} f s s+N S f_{s}-(M / K) f-(1-K) M f_{K}=0\tag{3} S2fss+NSfs(M/K)f(1K)MfK=0(3)
    敲黑板!!!BAW方法为何是半解析解,就在于我们要对式(3)进行简化,把最后一项省略掉,式(3)就可以变成一个二次常微分方程。那么最后一项为何可以省略掉呢?当商品期权距离到期日时间非常短(T=0)的时候, f K f_{K} fK接近于0,即执行价格对风险溢价几乎没影响。当商品期权距离到期日的时间非常长( τ = ∞ \tau=\infin τ=)的时候,K=1则1-K=0。因而, f f f可以通过近似求解以下的表达式得到:
    S 2 f s s + N S f s − ( M / K ) f = 0 S^{2} f s s+N S f_{s}-(M / K) f=0 S2fss+NSfs(M/K)f=0
    假设 f = a S q f=a S^{q} f=aSq,代入上式,可以得到两个特征根(M/K为正):
    q 1 = − ( N − 1 ) − ( N − 1 ) 2 + 4 M / K 2 < 0 q _ {1} = - \left( N - 1 \right) - \frac{\sqrt { \left( N - 1 \right) ^ { 2 } + 4 M / K }}{2}<0 q1=(N1)2(N1)2+4M/K <0

    q 2 = − ( N − 1 ) + ( N − 1 ) 2 + 4 M / K 2 > 0 q _ {2} = - \left( N - 1 \right) + \frac{\sqrt { \left( N - 1 \right) ^ { 2 } + 4 M / K }}{2}>0 q2=(N1)+2(N1)2+4M/K >0

    则有通解:
    f ( S ) = a 1 S q 1 + a 2 S q 2 , f \left( S \right) = a _ { 1 } S ^{ q _ { 1 } } + a _ { 2 } S ^{ q_2 } , f(S)=a1Sq1+a2Sq2,
    现在的问题是如何确定参数 a 1 a_1 a1 a 2 a_2 a2。这个问题要追溯到美式看涨期权的边值条件
    ∂ C ∂ S ( S ∗ , τ ) = 1 (B1) \frac{\partial C}{\partial S}(S^*,\tau)=1\tag{B1} SC(S,τ)=1(B1)

    C ( 0 , τ ) = 0 (B2) C(0,\tau)=0\tag{B2} C(0,τ)=0(B2)

    C ( S ∗ , τ ) = S ∗ − K (B3) C(S^*,\tau)=S^*-K\tag{B3} C(S,τ)=SK(B3)

    由(B2): C ( 0 , τ ) = c ( 0 , τ ) + K ( τ ) f ( 0 ) = 0 C(0,\tau)=c(0,\tau)+K(\tau)f(0)=0 C(0,τ)=c(0,τ)+K(τ)f(0)=0, 可以推出 f ( 0 ) = 0 f(0)=0 f(0)=0,因而 f ( S ) = a 2 S q 2 ( q 1 < 0 ) f(S)=a _ { 2 } S ^{ q_2 }(q_1<0) f(S)=a2Sq2(q1<0)

    由(B1): ∂ C ∂ S ( S ∗ , τ ) = ∂ c ∂ S ( S ∗ , τ ) + K ( τ ) ∂ f ∂ S \frac{\partial C}{\partial S}(S^*,\tau)=\frac{\partial c}{\partial S}(S^*,\tau)+K(\tau)\frac{\partial f}{\partial S} SC(S,τ)=Sc(S,τ)+K(τ)Sf 可以推导出
    1 = e ( b − r ) τ   N [ d 1 ( S ∗ ) ] + K q 2 a 2 S ∗ q 2 − 1 1=e^{(b-r) \tau} \mathrm{~N}\left[d_{1}\left(S^{*}\right)\right]+K q_{2} a_{2} S^{* q_{2}-1} 1=e(br)τ N[d1(S)]+Kq2a2Sq21
    这里用到了欧式期权的BS公式:
    c ( S , T ) = S e ( b − r ) T   N ( d 1 ) − X e − r T   N ( d 2 ) c(S, T)=S e^{(b-r) T} \mathrm{~N}\left(d_{1}\right)-X e^{-r T} \mathrm{~N}\left(d_{2}\right) c(S,T)=Se(br)T N(d1)XerT N(d2)
    其中 d 1 = [ ln ⁡ ( S / X ) + ( b + 0.5 σ 2 ) T ] / σ T , d 2 = d 1 − σ T d_{1}=\left[\ln (S / X)+\left(b+0.5 \sigma^{2}\right) T\right] / \sigma \sqrt{T}, d_{2}=d_{1}-\sigma \sqrt{T} d1=[ln(S/X)+(b+0.5σ2)T]/σT ,d2=d1σT

    由此,可以求解出
    a 2 = { 1 − e ( b − r ) τ   N [ d 1 ( S ∗ ) ] } / K q 2 S ∗ q 2 − 1 a_{2}=\left\{1-e^{(b-r) \tau} \mathrm{~N}\left[d_{1}\left(S^{*}\right)\right]\right\} / K q_{2} S^{* q_{2}-1} a2={1e(br)τ N[d1(S)]}/Kq2Sq21
    由(B3)可以求解出最优执行边界满足以下的表达式
    S ∗ − X = c ( S ∗ , τ ) + { 1 − e ( b − r ) τ   N [ d 1 ( S ∗ ) ] } S ∗ / q 2 S^{*}-X=c\left(S^{*}, \tau\right)+\left\{1-e^{(b-r) \tau} \mathrm{~N}\left[d_{1}\left(S^{*}\right)\right]\right\} S^{*} / q_{2} SX=c(S,τ)+{1e(br)τ N[d1(S)]}S/q2
    综上所述,美式看涨期权的解析公式如下所示:
    C ( S , τ ) = c ( S , τ ) + A 2 ( S / S ∗ ) q 2 ,  when  S < S ∗ C ( S , τ ) = S − K ,  when  S ≥ S ∗ \begin{array}{ll} C(S, \tau)=c(S, \tau)+A_{2}\left(S / S^{*}\right)^{q_{2}}, & \text { when } S<S^{*} \\ C(S, \tau)=S-K, & \text { when } S \geq S^{*} \end{array} C(S,τ)=c(S,τ)+A2(S/S)q2,C(S,τ)=SK, when S<S when SS
    其中 A 2 = ( S ∗ / q 2 ) { 1 − e ( b − r ) T   N [ d 1 ( S ∗ ) ] } A_{2}=\left(S^{*} / q_{2}\right)\left\{1-e^{(b-r) T} \mathrm{~N}\left[d_{1}\left(S^{*}\right)\right]\right\} A2=(S/q2){1e(br)T N[d1(S)]}

    美式看跌期权与看涨期权的求解套路差不多:
    ϵ p ( S , T ) = P ( S , T ) − p ( S , T ) \epsilon_p\left( S , T \right) = P \left( S , T \right) -p \left( S , T \right) ϵp(S,T)=P(S,T)p(S,T)
    区别在于边界条件不同
    ∂ P ∂ S ( S ∗ , τ ) = − 1 (B1) \frac{\partial P}{\partial S}(S^*,\tau)=-1\tag{B1} SP(S,τ)=1(B1)

    P ( ∞ , τ ) = 0 (B2) P(\infin,\tau)=0\tag{B2} P(,τ)=0(B2)

    P ( S ∗ , τ ) = K − S ∗ (B3) P(S^*,\tau)=K-S^*\tag{B3} P(S,τ)=KS(B3)

    python程序实现:参考Python implementation of the Barone-Adesi And Whaley model for the valuation of American options and their greeks. · GitHub

    调用实例:

    from BAW import *
    import matplotlib.pyplot as plt
    
    Stock = []
    Put_American = []
    Put_European = []
    Expir = []
    K = 50
    
    for S in range(1,151,10):
        Stock.append(S)
        Expir.append(max(K-S,0))
        put_American = getValue('American', 'Value', 'Put', S, K, 5, 0.03, 0.02,0.5)
        put_European = getValue('European', 'Value', 'Put', S, K, 5, 0.03, 0.02,0.5)
        Put_American.append(put_American)
        Put_European.append(put_European)
        
    plt.plot(Stock, Put_American,linewidth=0.5,label='American Option Value')
    plt.plot(Stock, Put_European,'r-.',linewidth=0.5,label='European Option Value')
    plt.plot(Stock, Expir, 'g--',linewidth=0.5)
    plt.legend()
    plt.show()
    

    [外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-5Y1Qch1q-1618749362203)(C:\Users\Cathy&Allen\Documents\草稿\BAW.png)]

    展开全文
  • 根据纳比吉安电磁法书上公式自己编写的长导线均匀空间解析解的C/C++代码,跟数值解做了对比,是吻合的,希望对你们有所帮助!
  • 将弹性空间地基受任意横向荷载作用下的静力位移积分变换解与两端自由梁的弯曲解析解相结合,采用三角级数展开的方法,对地基反力不做任何假设,求得了弹性空间地基上两端自由梁受任意横向荷载作用下的解析解,...
  • 提出一种计算变截面梁横向振动特性的半解析法。基于欧拉-伯努利梁理论给出的弯曲刚度、质量分布沿梁轴线连续或非连续变化的变截面梁横向振动方程;...通过与有限元法的数值结果比较说明半解析解的高精度和有效性。
  • 用波函数展开法给出了自由场地,将地球上一些峡谷用球形来模拟,推导出球形峡谷在平面SV波入射下的地表位移,总结出有利于实际工程应用的结论。空间表面的各点包括峡谷的边缘点的位移幅值没有表现出很高的...
  • 利用复变函数法,分析了部分裂纹面上受反平面剪应力和面内电载荷共同作用下有限高狭长压电体中含共线双无限裂纹问题,导出了电不可通边界条件下两个裂纹尖端场强度因子和机械应变能释放率的解析解.当不考虑电场作用...
  • 与基于弹性空间地基受任意竖向荷载作用的位移积分解建立的板与地基变形协调方程相结合,用三角级数法,得出弹性空间地基上四边自由正交异性矩形中厚板受任意竖向荷载作用的解析解。用该方法对算例进行计算,并将...
  • 采用双重Fourier变换,分析得到弹性空间地基受竖向稳态荷载作用下的积分变换解与四边自由矩形板的振动解析解相结合,得出弹性空间地基上四边自由矩形板稳态振动的解析解还给出算例及参数影响分析
  • AnaFlow为地下水流量方程式提供了几种分析和分析解决方案。 安装 您可以使用以下命令安装最新版本: pip install anaflow AnaFlow的文档 您可以在下找到文档。 例子 在下文中,众所周知的Theis函数被称为三个...
  • 认为桩体为渗透材料,给出混凝土管桩沉桩后桩周土中超静孔隙水压力消散的解析解,分析了孔压消散的因素影响。结果表明,土的黏弹性对孔压消散有明显影响,土的黏系数越大,超静孔压消散越快;在土的渗透性较小时桩身渗透...
  • 并与三个广义位移变量描述的弹性地基上四边自由正交各向异性矩形中厚板的弯曲控制方程相结合,用三角级数法,得出横观各向同性弹性空间地基上四边自由正交异性矩形中厚板受任意竖向荷载作用的弯曲解析解。...
  • 引入复位函数,对广义多极技术、多极理论和新型等效源方法等几种不同的半解析方法进行了统一描述,简化了半解析解的表述方式,降低了计算量,而且为直接从边值问题出发寻找复位函数解答建立了一条新途径,并进行了...
  • GD法...针对矩形板的动力响应问题,在空间域采用离散的GD法,在时间域取解析函数,构造了求解结构动力响应的GD半解析法。经实例计算表明,该方法是一种精度好效率高的求解动力响应问题的计算方法。
  • 结合虚拟激励法与Galerkin法,研究弹性圆拱在水平随机地震作用下随机响应的半解析解.在建立圆拱平面内动力平衡微分方程的基础上,通过选取适当的试函数,应用Galerkin法将动力平衡微分方程转化为线性常微分方程组....
  • GD法(General differential method)是从泰勒展开式...对新的控制方程在时间域取解析函数,在空间域采用离散的GD法,从而构造了卷积型GD半解析法。经对瞬态热传导问题的计算表明,该方法是一种精度好效率高的求解动力响应
  • 卷积型的Gurtin变分原理是目前在数学上唯一能和动力学初值问题完全等价的变分原理,它完全...对新的控制方程在时间域取解析函数,在空间域采用离散的GD法,从而构造了卷积型GD半解析法。该方法既可以达到和Gurtin变分原理
  • 永久美式期权解析解

    2021-04-30 21:31:59
    基本思想:有限行权期限的美式期权是没有解析解的,最多能找到一个半解析解,详见:美式期权的BAW定价方法 但是,如果我们假设无限到期日的美式期权(T→inf⁡T\to\infT→inf),即永久美式期权,则我们可以找到一

    永久美式期权解析解

    参考文章:Merton, R. C. . (1973). The theory of rational option pricing. The Bell Journal of Economics and Management Science, 4(1), 141-183.

    基本思想:有限行权期限的美式期权是没有解析解的,最多能找到一个半解析解,详见:美式期权的BAW定价方法

    但是,如果我们假设无限到期日的美式期权( T → inf ⁡ T\to\inf Tinf),即永久美式期权,则我们可以找到一个解析解。以美式看跌期权为例,假设 C C C是最优执行边界。则永久美式看跌期权的价格 P P P满足如下的表达式:
    1 2 σ 2 S 2 P s s + r S P s − r P = 0 , if C ≦ S ≦ ∞ (1) \frac{1}{2} \sigma^{2} S^{2} P_{ss}+r S P_{s}-r P=0,\quad\text{if}\quad C \leqq S \leqq \infty\tag{1} 21σ2S2Pss+rSPsrP=0,ifCS(1)
    subject to the boundary conditions
    P ( ∞ , ∞ ; K ) = 0 P ( C , ∞ ; K ) = K − C P s ( C , ∞ ; K ) = − 1 \begin{array}{l} P(\infty, \infty ; K)=0 \\ P(C, \infty ; K)=K-C\\ P_s(C, \infty ; K)=-1 \end{array} P(,;K)=0P(C,;K)=KCPs(C,;K)=1
    假设 K K K表示行权价格。与美式期权定价不同,由于这里的距离到期日的时间 τ = T − t \tau=T-t τ=Tt趋向于无穷到,(1)中并没有包含时间的偏导数项,因而公式(1)并非一个椭圆的偏微分方程,而是一个二次的常微分方程。它的通解可以表示成如下的形式:
    P ( S , ∞ ; K ) = a 1 S + a 2 S − γ (2) P \left( S,\infin;K \right) = a _ { 1 } S + a _ { 2 } S ^{ -\gamma }\tag{2} P(S,;K)=a1S+a2Sγ(2)
    其中, q 1 q_1 q1 q 2 q_2 q2为常微分方程(1)的特征解:


    假设 P = a S q P=aS^{q} P=aSq,代入到式(1)中,可得
    q 2 + ( γ − 1 ) q − γ = 0 q^2+(\gamma-1)q-\gamma=0 q2+(γ1)qγ=0
    其中 γ = 2 r / σ 2 \gamma=2r/\sigma^2 γ=2r/σ2。求解上式可得
    q 1 = 1 > 0 q 2 = − γ < 0 q_1=1>0\\ q_2=-\gamma<0 q1=1>0q2=γ<0


    由边界条件 P ( ∞ , ∞ ; K ) = 0 P(\infty, \infty ; K)=0 P(,;K)=0, 可知 a 1 = 0 a_1=0 a1=0。由边界条件 P ( C , ∞ ; K ) = K − C P(C, \infty ; K)=K-C P(C,;K)=KC, 可以求解出 a 2 a_2 a2 a 2 = ( K − C ) C γ a_2=(K-C)C^{\gamma} a2=(KC)Cγ

    因而我们有
    P ( S , ∞ ; K ) = ( K − C ) ( S / C ) − γ (2) P \left( S,\infin;K \right) = (K-C)(S/C)^{-\gamma} \tag{2} P(S,;K)=(KC)(S/C)γ(2)
    下面,我们来讨论如何确定最优执行边界 C C C。根据第三个边值条件| P s ( C , ∞ ; K ) = − 1 P_s(C, \infty ; K)=-1 Ps(C,;K)=1,可得
    C = γ K / ( 1 + γ ) C=\gamma K /(1+\gamma) C=γK/(1+γ)

    展开全文
  • 运用保角变换方法寻求了一个无限渗流域中非对称井流问题的解析解,分析了非对称井流水头势的分布规律.在此基础上,又利用井群工作时降深的叠加原理,进一步研究了非对称井群协同工作情况,获得了非对称井群共同作用下...
  • 然后根据模态叠加的方法构造悬臂梁的挠度函数,再利用虚功原理推导出集中力突然撤去情况下的自由阻尼悬臂梁瞬态响应近似解析解。算例分析表明:推导的公式准确可靠,且该方法简单,便于应用于工程计算。
  • 使用半解析环形板单元求解零部件转动时的固有频率(零部件转动时的固有频率称 为动频),推导出离心应力使单元增加的附加刚度矩阵,给出动频的求解过程,计算的结果被实验 证明。这种单元适用于研究金属切割锯片、盘形锥...
  • 10) %把轨道的周期和全周期标在图上 text(0.8*Ye(3,1),-2e7+Ye(3,2),['t3=' sprintf('%6.0f',Te(3))]) text(0.8*Ye(2,1),1.5e7+Ye(2,2),['t2=' sprintf('%6.0f',Te(2))]) %在x-y坐标上画地球 [XE,YE,ZE] = sphere...

    采用最简洁格式的ODE文件和解算指令,研究围绕地球旋转的卫星轨道。

    图 5.14.2.1-1 地球轨道卫星运动加速度

    (1)问题的形成

    轨道上运动的卫星,在Newton第二定律 ,和万有引力定律 作用下,有 。即 , ,而 。这里 是引力常数,

    是地球的质量。又假定卫星以初速度 在 处进入轨道。

    (2)构成一阶微分方程组

    令 ,则

    (5.14.2.1-5)

    初始条件为

    (5.14.2.1-6)

    (3)根据式(5.14.2.1-5)编写最简洁格式的ODE文件

    [dYdt.m]

    function Yd=DYdt(t,Y)

    % t 一定是标量形式的自变量

    % Y 必须是列向量

    global G

    ME %在函数中定义全局变量传递参数

    xy=Y(1:2);Vxy=Y(3:4); %前两个元素是"位移量",后两个是"速度量"。

    r=sqrt(sum(xy.^2));

    Yd=[Vxy;-G*ME*xy/r^3]; %Yd必须按式(5.14.2.1-5)编写,是与Y同维的列向量。

    (4)对微分方程进行解算

    global G

    ME %在主程序中定义全局变量传递参数 <1>

    G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;

    tspan=[t0,tf]; %指定解算微分方程的时间区间

    Y0=[x0;0;0;vy0]; %按式(5.14.2.1-6)给定初值向量

    [t,YY]=ode45('DYDt',tspan,Y0); %<8>

    X=YY(:,1); %输出Y的第一列是位移数据x(t)

    Y=YY(:,2); %输出Y的第二列是位移数据y(t)

    plot(X,Y,'b','Linewidth',2); hold on

    axis('image') %保证x、y轴等长刻度,且坐标框恰包容图形

    [XE,YE,ZE] =

    sphere(10); %产生单位球面数据(三维坐标)

    RE=0.64e7; %地球半径

    XE=RE*XE;YE=RE*YE;ZE=0*ZE; %坐标纸上的地球平面数据

    mesh(XE,YE,ZE),hold

    off %绘地球示意图

    图 5.14.2.1-2

    卫星轨道

    上例中,程序间的参数(如G和ME)传送,是依靠全局变量形式实现的。一般说来,编写程序时,应尽量少用全局变量,以免引起混乱。本例演示参数如何在指令间直接传送。

    要实现参数直接传送,必须对上例中的程序进行修改,具体如下:

    (1)重写ODE文件

    [DYDt2.m]

    function Yd=DYDt2(t,Y,flag,G,ME)

    % flag 按ODE文件格式规定,必须是第三输入宗量。对它的赋值由ode45指令自动产生。

    % 第4、5宗量是被传递的参数

    switch flag

    case

    '' %按规定:这里必须是空串。在此为"真"时,完成以下导数计算。

    X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V;-G*ME*X/r^3];

    otherwise

    end

    (2)按以下方法修改上例第(4)步中的程序,并运行之,可得相同结果。

    删去原主程序第<1>条指令;把原主程序的第<8>条指令改写为

    [t,YY]=ode45('DYDt2',tspan,Y0,[],G,ME); %第4宗量取缺省设置

    %第5、6宗量是被传递参数

    解算指令较复杂格式的使用示例

    带事件设置的ODE文件及主程序编写演示。本例将以较高精度计算卫星经过近地点和远地点的时间,并在图上标志。

    (1)ODE文件的编写

    [DYDt3.m]

    function varargout=DYDt3(t,Y,flag,G,ME,tspan,Y0)

    % DYDt3.m 供主程序调用的ODE函数文件

    % 本文件自带三个子函数:f,fi,fev。

    % t,

    Y 分别是自变量和一阶函数向量,是最基本的输入宗量。

    % flag 第三输入宗量,它专供解算指令(如ode45)作调用通知。

    % 在运行中,解算指令会根据需要向flag发不同的字符串。

    % varargout 是"变"输出宗量。它由变维的元胞数组构成。每个元胞中可以存放指令

    % 所产生的任意形式的数据。

    switch flag

    case

    '' % 必须用空串符。情况为"真"时,计算导数 dY/dt = f(t,Y)。

    varargout{1} =

    f(t,Y,G,ME); %输出为一个元胞,容纳f子函数的一个输出Yd

    case 'init' %

    必须用'init',情况为"真"时,传送计算区间、初值、设置参数。

    [varargout{1:3}] = fi(tspan,Y0);

    %输出为三个元胞,容纳fi子函数的三个输出量

    case 'events' % 必须用'events',情况为"真"时,设置事件性质。

    [varargout{1:3}] = fev(t,Y,Y0);

    %输出三个元胞,容纳fev子函数的三个输出量

    otherwise

    error(['Unknown flag ''' flag '''.']);

    end

    %

    ------------------------------------------------------------------

    function Yd = f(t,Y,G,ME)

    %计算导数子函数,被"父"函数DYDt3调用。

    X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V; -G*ME*X/r^3];

    %

    ------------------------------------------------------------------

    function [ts,y0,options] = fi(tspan,Y0)

    %设置时间区间、初值、算法参数子函数,被"父"函数DYDt3调用。

    ts=tspan;y0 = Y0;

    options = odeset('Events','on','Reltol',1e-5,'Abstol',1e-4);

    %开动ode45的"事件判断"功能,设置相对误差和绝对误差。

    %

    ------------------------------------------------------------------

    function [value,isterminal,direction] = fev(t,Y,Y0)

    %事件判断子函数,被"父"函数DYDt3调用。

    dDSQdt = 2 * ((Y(1:2)-Y0(1:2))' * Y(3:4));

    �SQdt是离初始点的距离u的时间导数du/dt,即u=(x-x0)^2+(y-y0)^2 。

    value = [dDSQdt; dDSQdt]; %定义两个穿越0的事件

    direction = [1; -1]; %第一事件:以渐增方式穿越0。第二事件:以渐减方式穿越0

    isterminal = [1; 0]; %第一事件发生后,终止计算;而第二事件发生后,继续计算。

    (2)运行以下主程序

    tspan=[t0,tf];Y0=[x0;0;0;vy0];

    [t,YY,Te,Ye,Ie]=ode45('DYDt3',[],[],[],G,ME,tspan,Y0); % <3>

    X=YY(:,1);Y=YY(:,2);

    plot(X,Y,'b','Linewidth',2);hold on

    text(0,6e7,'轨道','Color','b') %产生蓝色文字注释

    axis('image'); %保证x、y轴等长刻度,且坐标框恰包容图形

    %在三个事件发生点上画标记

    plot(Ye(1,1),0.4e7+Ye(1,2),'r^','MarkerSize',10)

    plot(Ye(2,1),0.4e7+Ye(2,2),'bv','MarkerSize',10)

    plot(Ye(3,1),-0.4e7+Ye(3,2),'b^','MarkerSize',10)

    %把轨道的半周期和全周期标在图上

    text(0.8*Ye(3,1),-2e7+Ye(3,2),['t3='

    sprintf('%6.0f',Te(3))])

    text(0.8*Ye(2,1),1.5e7+Ye(2,2),['t2='

    sprintf('%6.0f',Te(2))])

    %在x-y坐标上画地球

    [XE,YE,ZE] = sphere(10);RE=0.64e7;XE=RE*XE;YE=RE*YE;ZE=0*ZE;

    mesh(XE,YE,ZE)

    text(1e7,1e7,'地球','Color','r'), hold

    off %产生红色文字注释

    图 5.14.2.2-1 带事件标注的卫星轨道图

    关于ODE文件的说明

    (1)ODE文件的模板

    下面就是MATLAB的模板文件odefile.m ,为便于读者阅读,本书作者适当地给以注解说明。

    [odefile.m]

    function varargout = odefile(t,y,flag,p1,p2)

    % odefile.m 是MATLAB提供的模板文件。用户可根据需要,使用任何其他文件名。

    % varargout 是MATLAB专门设计的"变长度"输出宗量。它由元胞数组构成,因此

    % 可适应任意多个、任意形式的输出宗量。请参看第8.5.2节。

    %

    t 自变量。不管在下面指令中是否出现,自变量t必须作为第一输入宗量。

    %

    y 一阶微分方程组的列向量形式函数。它必须作为第二输入宗量。

    %

    flag 切换变量。它必须处在第三输入宗量位置。用户无须也不要对它直接赋

    % 值;它的赋值由微分方程解算指令(solver),自动产生。

    % p1,

    p2 是被传递的参数。这里,作为示意,仅列出两个。

    % 根据需要,传递参数的数目不受限制。

    % 以下是switch-case构成的多向选择控制。应注意:

    % (A)各情况的情况表达式是MATLAB规定的,不得任意改变。

    % (B)各情况所执行的任务也是规定的,不得任意改变。

    % (C)各情况内,变长输出宗量元胞数是规定的,不得任意改变。

    % (D)各情况内,(等式右边的)函数名称可以改变,但相应子函数名称要一致。

    % (E)各情况内,(等式右边的)函数中所包含的t , y 是必须的,

    % 不管它们是否以显式出现相应子函数体中。

    % (F)各情况内,(等式右边的)函数中所包含的p1, p2

    是传递参数的示意性表示。

    % 具体参数视需要而定,只要该参数已经传入odefile.m 函数内存。

    switch flag

    case

    '' % 规定空字符串 ''

    情况:专管一阶导数 dy/dt = f(t,y)的计算

    varargout{1} =

    f(t,y,p1,p2);

    case 'init' % 规定 'init' 情况:专管三个宗量 [tspan ,

    y0 , options]的设置。

    [varargout{1:3}] =

    init(p1,p2);

    case 'jacobian'

    % 规定

    'jacobian' 情况:专管计算解析的 Jacobian 矩阵 df/dy。

    varargout{1} =

    jacobian(t,y,p1,p2);

    case 'jpattern'

    % 规定

    'jpattern' 情况:专管计算稀疏的数值 Jacobian 矩阵 df/dy。

    varargout{1} =

    jpattern(t,y,p1,p2);

    case

    'mass '

    % 规定 'mass'

    情况:专管计算质量矩阵(mass matrix)。

    varargout{1} =

    mass(t,y,p1,p2);

    case 'events' %

    规定'events'情况:专管事件定义和判断

    [varargout{1:3}] =

    events(t,y,p1,p2);

    otherwise

    error(['Unknown flag '''

    flag '''.']);

    end

    % 以下是odefile.m 的子函数,它们分别与"父"函数中各情况对应。

    %---------------------------------------------------------------------------

    function dydt = f(t,y,p1,p2)

    % 空字符串 '' 情况内调用的子函数:专管一阶导数 dy/dt的计算

    % dydt 一阶导数。该变量名可由用户按需要取名。

    % 下面函数内容,由用户自己编写,最后产生输出dydt 。

    dydt =〈由用户编写〉

    %

    ---------------------------------------------------------------------------

    function [tspan,y0,options] =

    init(p1,p2)

    % 'init' 情况内调用的子函数:专管三个宗量 [tspan , y0 , options]的设置。

    % 三个输出宗量的名称可由用户自起,但各位置上宗量的性质不能改变。

    % 下面函数内容,由用户自己编写,最后产生输出三个输出 。

    tspan =〈积分的时间区段或时间点集〉

    y0 =〈初值〉

    options =〈或用 odeset 设置算法参数,或用空阵符 [ ]

    表示缺省设置〉

    %

    ---------------------------------------------------------------------------

    function dfdy =

    jacobian(t,y,p1,p2)

    % 'jacobian' 情况调用的子函数:专管计算解析的 Jacobian 矩阵 df/dy。

    % 输出宗量的名称可由用户自起。

    % 下面函数内容,由用户自己编写,最后产生输出输出dfdy 。

    dfdy =〈 Jacobian 矩阵〉

    %

    ---------------------------------------------------------------------------

    function S = jpattern(t,y,p1,p2)

    % 'jpattern' 情况调用的子函数:专管计算解析的 Jacobian 矩阵 df/dy。

    % 输出宗量的名称可由用户自起。

    % 下面函数内容,由用户自己编写,最后产生输出S 。

    S =〈稀疏形式的数值 Jacobian 矩阵〉

    %

    ---------------------------------------------------------------------------

    function M = mass(t,y,p1,p2)

    % 'mass' 情况调用的子函数:专管计算质量矩阵(mass matrix)。解刚性方程时用。

    % 输出宗量的名称可由用户自起。

    % 下面函数内容,由用户自己编写,最后产生输出M 。

    M =〈质量矩阵〉

    %

    --------------------------------------------------------------------------

    function [value,isterminal,direction] =

    events(t,y,p1,p2)

    % 'events' 情况内调用的子函数:专管三个宗量 [value,isterminal,direction]的设置。

    % 三个输出宗量的名称可由用户自起,但各位置上宗量的性质不能改变。

    % 下面函数内容,由用户自己编写,最后产生输出三个输出 。

    value =〈各元素是标量事件函数,即函数穿越0点为事件〉

    direction =〈各元素定义事件函数穿越方式:1为向正穿越0;-1为向负穿越0;0不管方向〉

    isterminal =〈各元素定义事件发生后计算是否继续:1为终止计算;0为继续计算〉

    %

    --------------------------------------------------------------------------

    关于解算指令选项options的属性设置

    options属性处理和输出函数使用演示

    仍以卫星轨道问题为例。本例演示如何通过对options域的直接设置,借助微分方程解算输出指令,表现解算的中间结果。具体目标是:画出解向量

    中由 构成的相平面。相平面的绘制是在微分方程解算中间逐步完成的。

    (1)编写ODE文件DYDt4.m (在DYDt3.m 基础上删改而成)

    [DYDt4.m]

    function varargout=DYDt4(t,Y,flag,G,ME,tspan,Y0)

    switch flag

    case ''

    varargout{1} =

    f(t,Y,G,ME);

    case 'init'

    [varargout{1:3}] = fi(tspan,Y0);

    otherwise

    error(['Unknown flag ''' flag '''.']);

    end

    %

    ------------------------------------------------------------------

    function Yd = f(t,Y,G,ME)

    X=Y(1:2);V=Y(3:4);r=sqrt(sum(X.^2));Yd=[V; -G*ME*X/r^3];

    %

    ------------------------------------------------------------------

    function [ts,y0,options] = fi(tspan,Y0)

    ts=tspan;y0 = Y0;

    % 采用向域直接赋值法,设置options属性。以供与odeset使用方法对照。

    options.RelTol=1e-5;options.AbsTol=1e-4;

    options.OutputFcn='odephas2'; %在积分进程中,绘制相平面图。

    options.OutputSel=[1 3];%解向量的第1、3分量分别为相平面图的横、纵坐标量。

    (2)编写脚本文件odeexp4.m

    [odeexp4.m]

    %odeexp4.m

    G=6.672e-11;ME=5.97e24;vy0=4000;x0=-4.2e7;t0=0;tf=60*60*24*9;

    tspan=[t0,tf];Y0=[x0;0;0;vy0];

    % 以下四条指令,为相平面图预置一个坐标轴范围

    clf,set(gca,'xlim',[-5 25]*1e7,'ylim',[-3

    3]*1e3);%设置适当坐标范围 <5>

    box

    on %框形坐标

    hold

    on; %使中间结果绘在同一幅图上

    ode45('DYDt4',[],[],[],G,ME,tspan,Y0);hold

    off

    (3)在MATLAB指令窗中运行以下指令,即可在图形窗中见到相平面图的逐步绘制。

    shg,odeexp4

    图 5.14.4.2-1 微分方程输出函数相平面图的实时绘制

    展开全文
  • 采用假定的大圆弧面来模拟空间,利用波函数展开法和有限项Fourier级数展开技术,在频域内给出圆弧形沉积谷地对Rayleigh波三维散射的一个解析解。以此解答为基础,进一步分析入射波的频率和入射角度以及谷地深度对圆弧...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 24,063
精华内容 9,625
关键字:

半解析解