三次样条插值函数的Matlab代码

2019-01-03 16:24

并针对下面一组具体实验数据 0.25 0.5000 0.3 0.5477 0.39 0.6245 0.45 0.6708 .

0.53 0.7280 求解,其中边界条件为

解:Matlab计算程序为:

clear clc

x=[0.25 0.3 0.39 0.45 0.53]

y=[0.5000 0.5477 0.6245 0.6708 0.7280] n=length(x); for i=1:n-1

h(i)=x(i+1)-x(i); end

for i=1:n-2

k(i)=h(i+1)/(h(i)+h(i+1)); u(i)=h(i)/(h(i)+h(i+1)); end

for i=1:n-2

gl(i)=3*(u(i)*(y(i+2)-y(i+1))/h(i+1)+k(i)*(y(i+1)-y(i))/h(i)); end

g0=3*(y(2)-y(1))/h(1);

g00=3*(y(n)-y(n-1))/h(n-1); g=[g0 gl g00]; g=transpose(g) k1=[k 1]; u1=[1 u];

Q=2*eye(5)+diag(u1,1)+diag(k1,-1) m=transpose(Q\\g) symsX; for i=1:n-1

p1(i)=(1+2*(X-x(i))/h(i))*((X-x(i+1))/h(i))^2*y(i); p2(i)=(1-2*(X-x(i+1))/h(i))*((X-x(i))/h(i))^2*y(i+1); p3(i)=(X-x(i))*((X-x(i+1))/h(i))^2*m(i); p4(i)=(X-x(i+1))*((X-x(i))/h(i))^2*m(i+1); p(i)=p1(i)+p2(i)+p3(i)+p4(i); p(i)=simple(p(i)); end s1=p(1) s2=p(2) s3=p(3) s4=p(4) for k=1:4

for z=x(k):0.001:x(k+1)

q=eval(subs(p(k),'X','z')); plot(z,q,'b') holdon end end gridon legend('èy′??ùì??ú??') title('2??μ') xlabel('x') ylabel('p')

计算结果为:

s1 =

-(705394867539368680*X^3 - 529046150654530286*X^2 + 23087199381998953*X - 40023205577025431)/112589990684262400 s2 =

(257361898089296225*X^3)/136796838681378816 - (160081743506404901*X^2)/60798594969501696 + (404209705972252727*X)/202661983231672320 + 429142243010323951/3166593487994880000 s3 =

- (3495912536773825*X^3)/7599824371187712 + (1437374409830143*X^2)/13510798882111488 + (10427488839800859*X)/11258999068426240 + 12358231431982259943/45035996273704960000 s4 =

(38626753769033575*X^3)/18014398509481984 - (245666153971053021*X^2)/72057594037927936 + (3614707928905781673*X)/1441151880758558720 + 26732501105874704913/720575940379279360000

通过绘制曲线,发现为3次样条曲线,且数据拟合较好。


三次样条插值函数的Matlab代码.doc 将本文的Word文档下载到电脑 下载失败或者文档不完整,请联系客服人员解决!

下一篇:企业工会经费财务管理办法

相关阅读
本类排行
× 注册会员免费下载(下载后可以自由复制和排版)

马上注册会员

注:下载文档有可能“只有目录或者内容不全”等情况,请下载之前注意辨别,如果您已付费且无法下载或内容有问题,请联系我们协助你处理。
微信: QQ: