微分方程数值求解

数学星空

TO mathe:

$tu^2u''=-(1+u')^3+2tu(1+u')^2-(u+t)u(1+u')+tu^3$

做变量替换$v=\frac{1+u'}{u}$,得出微分方程$\frac{dv}{d\theta}=(\theta-v)(v^2+1)$ 似乎应改成

$\frac{dv}{d\theta}=(1-\frac{v}{\theta})(v^2+1)$

mathe

有道理。其实在v离1比较远时的图像是类似的，就在0附近区别比较大。而这个方程应该有v(0)=0,v'(0)=1/2,而且可以看成这个解是一个奇函数
设其在0的展开式为$v(t)=\sum_{n=0}^{+\infty} a_{2n+1}t^{2n+1}$,我们得出
$\sum_{n=0}^{+\infty}(2n+1)a_{2n+1}t^{2n+1}=(t-\sum_{n=0}^{+\infty} a_{2n+1}t^{2n+1})(1+(\sum_{n=0}^{+\infty} a_{2n+1}t^{2n+1})^2)$
展开比较系数就可以得出$a_1=1/2,(2n+2)a_{2n+1}=\sum_{k=0}^{n-1}a_{2k+1}a_{2n-2k-1}-\sum_{k=0}^{n-1}a_{2k+1}\sum_{m=0}^{n-k-1}a_{2m+1}a_{2(n-k-m)-1},n>1$
而对应的方程$u'-vu+1=0$初始条件$u_0(0)=0$的解也是奇函数。而$V(t)=exp(\int_0^t v(t)dt)$是偶函数，u的通解为$u(t)=u_0(t)+C*V(t)$

mathe

更新v得到结果如下 d = 0.59090489185668802875786709583348871021050937452181485289441819749025679734298579253579
a = 3.54925966919232488335307997505516088581017092471241443348333130503956135888576020826237
这个结果同fans的数据一致了
更新图片
(这是v关于$\theta$的函数)
(这是u关于$\theta$的函数)
(这是目标函数关于d的导数的图像)
(这是目标函数关于d的图像，也就是距离的期望值)
(这个是最优解对应的图像)
添加附件包含Pari/gp源代码
v,u函数在$\theta=0$的收敛半径好像大于$pi$,所以实际上，分段数目可以少很多，但是由于直接从以前代码修改过来，就保留了更多的分段了。
而比较有意思的是计算出来在$\theta=0.5$的收敛半径好像小于1，这个有点不符合道理（圆的边界上必然有函数的奇点，但是落在0为中心收敛圆内部）。
一个比较合理的解释就是计算误差引起的。可以看到满足v的微分方程的一族曲线中，除了这条v(0)=0的曲线，其余的在靠近$\theta=0$时，值都会发生剧烈变化，远远离开原点。
这个导致咋$\theta=0.5$时，如果有一点计算误差，对应到$\theta=0$处，会有很大的误差，所以就远离我们的目标曲线了。所幸，$\theta$越大，这种误差影响越小。
另外，如果我们能够修改一下代码，更大范围使用在$\theta=0$处的展开式(可以有精确表达形式),会有更好的精度。
另外比较有意思的是取$u_0(0)=0$对应的u的解在$\theta=0$处展开式的前面15项竟然全部是正数(对应到$\theta^31$系数)，但是可惜后面的就出现负数了。

wayne

七年了。
mathe竟然没有 七年之痒，

数学星空

算出u(x)后，d如何求得，然后a的表达式是什么呢？

mathe

看45#

wayne

为啥我算的跟mathe的不同：

1. ans=Series[v[t],{t,0,20}];
   
2. f=Normal[ans/.First@SolveAlways[Join[{t D[#,t]==(t-#)(#^2+1)&[ans]},{v[0]==0}],t]]

复制代码

$v[t] = t/2 + t^3/32 + t^5/768 - t^7/49152 - t^9/131072 - (
7 t^11)/14155776 + (19 t^13)/6341787648 + (
5179 t^15)/1623497637888 + (12715 t^17)/50096498540544 + (
3571 t^19)/1753377448919040$

数学星空

To wayne:

因为mathe在42#算的微分方程是

微分方程$\frac{dv}{d\theta}=(\theta-v)(v^2+1)$, 边界条件$v(0)=0$

条件下的解,当然不一样

后来mathe重新计算了,只是没有写出来..

数学星空

我还是没有弄明白mathe在43#提出的计算特解\(u_0(x)\)的方案:
下一步可以挑选一个合适的d,计算对应的u(x),由于我们已经有$u(2\pi-2d)=\tan(d)$,直接根据微分方程$u'(x)-v(x)u(x)+1=0$可以直接根据$v(x)$在$2\pi-2d$的展开式计算出$u(x)$在$2\pi-2d$的展开式，得到一个特解$u_0(x)$。

wayne

的确如mathe所言，分母含有大量的小素因子。v(t) 的49次幂 的系数 是$- (536894642831175368956914839 t^49)/634244296097068809167943077627551535953320974745600000$
分母分解起来是 {{2, 112}, {3, 18}, {5, 5}, {7, 3}, {11, 3}, {13, 1}, {17, 1}}
即 112 个2，18个3，5个5，3个7