应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜    

CyRhmU.jpeg
小木虫论坛-学术科研互动平台 » 专业学科区 » 数学 » 计算数学 » Laplace数值反演wooden方法请大家帮忙（如能解决，可加金币）
81/1返回列表
查看: 1622  |  回复: 7
只看楼主@他人 存档 新回复提醒 (忽略) 收藏    在APP中查看

皓小天

木虫之王 (职业作家)

[求助] 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不太熟悉。

希望在这方面擅长的虫友们帮我看看可不可以解决这个问题
拜谢Laplace数值反演wooden方法请大家帮忙（如能解决，可加金币）
2015-10-20_10-13-16.png
回复此楼

» 本帖附件资源列表

» 本帖@通知

@月只蓝

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:
请不要站内找我要书，如果需要请到书籍板块求助
1楼2015-10-20 10:33:10
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nagami

木虫 (正式写手)

【答案】应助回帖

...
V=zeros(1,N);
....

V(i)=V(i)+(-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);
...
这样写就ok了
要是大规模计算，这种写法效率差，应该避免不必要的重复性计算
赞一下 回复此楼
女靠衣装；男靠金装
2楼2015-10-23 13:51:21
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nagami

木虫 (正式写手)

【答案】应助回帖

...
V=zeros(1,N);
....

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));
...
是这样，0偏置和1偏置的区别，fortran和c/c++，起始数组编号不同
赞一下 回复此楼
女靠衣装；男靠金装
3楼2015-10-23 13:54:25
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

皓小天

木虫之王 (职业作家)
引用回帖:
3楼: Originally posted by nagami at 2015-10-23 13:54:25
...
V=zeros(1,N);
....

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));
...
是这样，0偏置和1偏置的区别，fortran和c/c++，起始数 ...

计算结果不对啊，你看我按照你的修改如下：
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明显不等啊
**********************************************
谢谢，先送你鲜花，在看看哪不对
赞一下(1人) 回复此楼
请不要站内找我要书，如果需要请到书籍板块求助
4楼2015-10-23 14:06:45
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

nagami

木虫 (正式写手)

【答案】应助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ...
皓小天: 金币+150, ★★★★★最佳答案, 非常感谢 2015-10-23 16:20:11
引用回帖:
4楼: Originally posted by 皓小天 at 2015-10-23 14:06:45
计算结果不对啊，你看我按照你的修改如下：
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)/(factor ...

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
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

358065752

新虫 (初入文坛)
你在VB中试过吗
回复此楼
6楼2017-03-20 20:49:42
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

遗忘消失的灆
7楼2017-03-31 16:21:09
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

camilali

新虫 (初入文坛)
楼主，得出结果了么？

发自小木虫Android客户端
赞一下 回复此楼
8楼2017-10-23 08:15:02
 已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 皓小天 的主题更新
81/1返回列表
普通表情 龙 兔 虎 猫
信息提示
关闭
请填处理意见
关闭