找回密码
 欢迎注册
楼主: mathe

[求助] 微分方程数值求解

[复制链接]
发表于 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)$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 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)$
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2015-1-22 08:16:17 来自手机 | 显示全部楼层
更新v得到结果如下 d = 0.59090489185668802875786709583348871021050937452181485289441819749025679734298579253579
a = 3.54925966919232488335307997505516088581017092471241443348333130503956135888576020826237
这个结果同fans的数据一致了
更新图片
v.gif (这是v关于$\theta$的函数)
u.gif (这是u关于$\theta$的函数)
dT-d.gif (这是目标函数关于d的导数的图像)
T-d.gif (这是目标函数关于d的图像,也就是距离的期望值)
l.gif (这个是最优解对应的图像)
添加附件包含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$系数),但是可惜后面的就出现负数了。

a.zip

2.19 KB, 下载次数: 2, 下载积分: 金币 -1 枚, 经验 1 点, 下载 1 次

评分

参与人数 2威望 +24 金币 +36 贡献 +24 经验 +24 鲜花 +24 收起 理由
KeyTo9_Fans + 12 + 20 + 12 + 12 + 12 神一般的存在
数学星空 + 12 + 16 + 12 + 12 + 12 神一般的存在

查看全部评分

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-1-22 12:34:54 | 显示全部楼层
七年了。
mathe竟然没有 七年之痒,

点评

代码贴出来了,貌似不太一样。 有空再细看  发表于 2015-1-22 22:50
用数学软件应该可以给出级数表示式,你试试看…  发表于 2015-1-22 13:37
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-1-22 13:47:34 | 显示全部楼层
算出u(x)后,d如何求得,然后a的表达式是什么呢?
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2015-1-22 14:12:31 来自手机 | 显示全部楼层
看45#
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
回复

使用道具 举报

发表于 2015-1-22 22:48:24 | 显示全部楼层
为啥我算的跟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$

点评

的确,看似很大,分解起来都是小因子,2的个数很多~  发表于 2015-1-23 12:57
这个没有错,我们可以发现分母具有大量小素因子  发表于 2015-1-23 06:19
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2015-1-22 23:41:26 | 显示全部楼层
To wayne:

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

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

条件下的解,当然不一样

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

点评

哦 哦 。我还以为 mathe 前面已经都更新了  发表于 2015-1-23 12:58
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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)$。

点评

另外在得出v在0的展开式以后,也可以直接根据微分方程给出u(0)=0的特解,这个的好处可以同样得到一个具有精确表达式的u(在0附近)  发表于 2015-1-23 06:19
你可以参考54#中我的Pari/Gp代码,函数gen_u_params中代码: e=2*Pi-2*d; a0=fun_v(v_params,e);a=fun(a0,e,n); b=gen_u_from_v(tan(d),a0,a);m0=use_fun(tan(d),e,b,5.0); 就是解u在这里的展开式。  发表于 2015-1-23 06:17
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 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

点评

令 t= v[t], 计算一下  发表于 2015-1-23 18:39
另外v函数在复平面半径3.4左右的圆上有奇点,而且应该是本性奇点。但是不知道怎么去计算  发表于 2015-1-23 17:23
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|数学研发网 ( 苏ICP备07505100号 )

GMT+8, 2024-3-28 21:09 , Processed in 0.047054 second(s), 18 queries .

Powered by Discuz! X3.5

© 2001-2024 Discuz! Team.

快速回复 返回顶部 返回列表