24小时热门版块排行榜

返回列表

huameitang05

银虫 (初入文坛)

应助: 0 (幼儿园)
金币: 250.7
帖子: 47
在线: 59.6小时
虫号: 99774

[交流] 【求助】急求帮忙解决MATLAB 样条差值求解微分代码出现的问题，使代码正常运行

function y0=hillsplines3(x0,x,y)
if size(x,1)~=1
x=x';
end
if size(y,1)~=1
y=y';
end
[a,t]=hillsplinex(x,y);a=a';t=t';
for i=1:length(t)-4
   phi(i)=powerplus3(x0,t(i:i+4));
end
y0=sum(a.*phi);
function [ax,tx]=hillsplinex(t,y)
h=diff(t);delta=diff(y)./h;n=length(h)+1
d1=pchipendpoint(h(1),h(2),delta(1),delta(2));
dn=pchipendpoint(h(n-1),h(n-2),delta(n-1),delta(n-2));
r=[d1;y';dn];
tx=[t(1)-[1:3]*h(1),t,t(length(t)+[1:3])*h(length(h))];
a=zeros(length(tx)-4,length(tx)-4);
for i=1:3
a(1,i)=depowerplus3(t(1),tx(i:i+4));
end
for i=length(tx)-6:length(tx)-4
a(length(tx)-4,i)=depowerplus3(t(length(t)),tx(i:i+4));
end
for j=2:length(tx)-5
   for i=j-1:j+1
      a(j,i)=powerplus3(t(j-1),tx(i:i+4));
   end
end
ax=a\r;tx=tx';
function y=powerplus3(x,t)
      c=t;
      for i=1:5
         c(i)=[];
         beta(i)=24/prod(t(i)-c);
         c=t;
      end
      powerplus=abs(x-t).^3;
      y=0.5*sum(beta.*powerplus);
      function y=depowerplus3(x,t)
         c=t;
      for i=1:5
         c(i)=[];
         beta(i)=24/prod(t(i)-c);
         c=t;
      end
      powerplus=3*sign(x-t).*(x-t).^2;
      y=0.5*sum(beta.*powerplus);
         function d=pchipendpoint(h1,h2,del1,del2)
         d=((2*h1+h2)*del1-h1*del2)/(h1+h2);
         if sign(d)~=sign(del1)
            d=0;
         elseif(sign(del1)~=sign(del2))&(abs(d)>abs(3*del1))
            d=3*del1;
         end

[ Last edited by huameitang05 on 2011-1-4 at 17:21 ]

回复此楼

» 猜你喜欢

» 抢金币啦！回帖就可以得到:

查看全部散金贴

1楼 2011-01-04 17:13:22

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

lijinfeng042

木虫 (小有名气)

仿真EPI: 2
应助: 1 (幼儿园)
金币: 2156.1
帖子: 291
在线: 31.5小时
虫号: 1019062

★ ★ ★ ★ ★
huameitang05(金币+1):谢谢参与
huameitang05(金币+1):谢谢意见 2011-01-04 19:14:27
xiegangmai(金币+2):辛苦了！ 2011-01-07 22:06:47
robert2020(金币+2):辛苦了！最近比较忙，评帖不及时，敬请见谅！ 2011-01-11 10:47:08

引用回帖:

Originally posted by huameitang05 at 2011-01-04 17:13:22:
function y0=hillsplines3(x0,x,y)
if size(x,1)~=1
x=x';
end
if size(y,1)~=1
y=y';
end
[a,t]=hillsplinex(x,y);a=a';t=t';
for i=1:length(t)-4
phi(i)=powerplus3(x0,t(i:i+4));
en ...

具体太多代码没看...只是说说思路离散数据的微分问题
既然你用的是样条插值那就说说这个方法
数据->csapi(x,y)->fnder(cs)->fnval(pp,x)