数学星空 发表于 2015-1-21 19:22:04

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 发表于 2015-1-21 21:39:10

有道理。其实在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 发表于 2015-1-22 08:16:17

更新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 发表于 2015-1-22 12:34:54

七年了。
mathe竟然没有 七年之痒, :lol

数学星空 发表于 2015-1-22 13:47:34

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

mathe 发表于 2015-1-22 14:12:31

看45#

wayne 发表于 2015-1-22 22:48:24

为啥我算的跟mathe的不同:
ans=Series,{t,0,20}];
f=Normal==(t-#)(#^2+1)&},{v==0}],t]]

$v = 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$

数学星空 发表于 2015-1-22 23:41:26

To wayne:

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

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

条件下的解,当然不一样

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

数学星空 发表于 2015-1-23 00:24:08

我还是没有弄明白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 发表于 2015-1-23 13:06:41

的确如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
页: 1 2 3 4 5 [6] 7 8
查看完整版本: 微分方程数值求解