| 查看: 1726 | 回复: 7 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
皓小天木虫之王 (职业作家)
|
[求助]
Laplace数值反演wooden方法请大家帮忙（如能解决，可加金币） 已有1人参与
|
||
|
希望能与各位喜欢数学的同道们讨论下。 wooden数值反演方法是stehfest方法的一种改进(常称为AWG方法)，在孔祥言的书中写出了其数值反演的具体方式(在下面图片中)，我带入matlab计算，发现反演的结果明显不对（代码在下面）。我将wooden原文献的代码（好像要fortran编写的）放在附件中，孔祥言应该是根据wooden的代码写的通式。孔祥言的通式如下（下面图片中）： 编写的matlab代码如下（按照孔祥言通式编写）： t=2; s=0; v=0; N=20; V=ones(1,N); f=@(x) 1./(1+x);%laplace空间的函数,很容易知道此函数的真实空间函数为exp(-t) for i=1:N for k=floor((i+1)/2):min(i,N/2) V(i)=v+(-1).^(N/2+i).*k.^(N/2)*factorial(2*k+1)/factorial(k+1)/factorial(k)/factorial(N/2-k+1)/factorial(i-k+1)/factorial(2*k-i+1); end s=s+log(2)./t.*V(i).*f(log(2)./t*i); end ******************************************* ******************************************* wooden的程序在附件中（wooden并没有写通式，在他文章中只有附件中的程序，还有与其他方法计算的对比结果，所以信息就只有这么多） 希望得到的东西：可以正确的运行wooden的程序，用MATLAB即可，得到正确的反演结果。 我现在不知道是孔祥言通式有问题，还是我编程出问题了（我认为是正确的），当然方法的来源都是来自wooden的程序，但是我对fortran不太熟悉。 希望在这方面擅长的虫友们帮我看看可不可以解决这个问题 拜谢   2015-10-20_10-13-16.png |
回复此楼» 本帖附件资源列表
-
欢迎监督和反馈：小木虫仅提供交流平台，不对该内容负责。
本内容由用户自主发布，如果其内容涉及到知识产权问题，其责任在于用户本人，如对版权有异议，请联系邮箱：xiaomuchong@tal.com - 附件 1 : 2015-10-20_10-13-16.png
- 附件 2 : Well_Test_Analysis_Benefits_From_New_Method_Of_Laplace_Space_Inversion.pdf
2015-10-20 10:13:39, 562.6 K
2015-10-20 10:20:40, 303.65 K
» 本帖@通知
@月只蓝
» 猜你喜欢
295材料工程专硕求调剂
已经有8人回复
-
求调剂
已经有8人回复
-
293求调剂
已经有5人回复
-
279求调剂
已经有8人回复
-
328求调剂
已经有8人回复
-
322求调剂：一志愿湖南大学 材料与化工（085600），已过六级。
已经有9人回复
-
307求调剂
已经有17人回复
-
286求调剂
已经有5人回复
-
085600 295分求调剂
已经有6人回复
-
287求调剂
已经有13人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- Matlab 7.0 基础教程 清华大学
已经有29人回复
【答案】应助回帖
★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
皓小天: 金币+150, ★★★★★最佳答案, 非常感谢 2015-10-23 16:20:11
皓小天: 金币+150, ★★★★★最佳答案, 非常感谢 2015-10-23 16:20:11
|
t=2; s=0; v=0; N=18; V=zeros(1,N); f=@(x) 1./(1+x);%laplace空间的函数,很容易知道此函数的真实空间函数为exp(-t) for i=1:N for k=floor((i+1)/2):min(i,N/2) V(i)=V(i) + k^(N/2)*factorial(2*k)/(factorial(k)*factorial(k-1)*factorial(N/2-k)*factorial(i-k)*factorial(2*k-i)); end V(i)=(-1)^(N/2+i)*V(i); s=s+log(2)./t.*V(i).*f(log(2)./t*i); end s exp(-t) |
5楼2015-10-23 15:43:23
2楼2015-10-23 13:51:21
3楼2015-10-23 13:54:25
皓小天
木虫之王 (职业作家)
- 应助: 19 (小学生)
- 贵宾: 0.1
- 金币: 76016.6
- 散金: 502
- 红花: 137
- 帖子: 3801
- 在线: 869小时
- 虫号: 1864745
- 注册: 2012-06-19
- 性别: GG
- 专业: 力学中的基本问题和方法
|
计算结果不对啊，你看我按照你的修改如下： t=2; s=0; v=0; N=20; V=zeros(1,N); f=@(x) 1./(1+x); for i=1:N for k=floor((i+1)/2):min(i,N/2) V(i)=V(i) + k^(N/2)*factorial(2*k)/(factorial(k)*factorial(k-1)*factorial(N/2-k)*factorial(i-k)*factorial(2*k-i)); end s=s+log(2)./t.*V(i).*f(log(2)./t*i); end ****************************** 结果如下： s=452409033762.281; **************************** laplace空间1/(1+x)对应的原函数为exp(-t),所以理论解在t=2的时候， exp(-2）=0.1353 这与452409033762.281明显不等啊 ********************************************** 谢谢，先送你鲜花，在看看哪不对 |
4楼2015-10-23 14:06:45
