計算方法大作業第二次

2021-03-04 09:41:17 字數 4684 閱讀 9568

數值計算第二次大作業

1.給定插值條件如下:

i 0 1 2 3 4 5 6 7

xi 8.125 8.4 9.0 9.485 9.6 9.959 10.166 10.2

yi 0.0774 0.099 0.280 0.60 0.708 1.200 1.800 2.177

作三次樣條函式插值,取第一類邊界條件

y0』=0.01087 y7』=100

根據題目要求,首先要構造三次樣條函式,三次樣條函式的構造過程如下:

設有<<…<共n個插值節點,任意給定一組常數,,…,,要求構造乙個插值三次樣條函式,使得如下插值條件得以滿足:

,i=0,1,…,n

經過插值點的三次樣條函式是一組三次多項式,即有:

由節點處的連續性可知:

由節點處的一階與二階光滑性可知:

又設,記,則。再根據邊界條件,從而可相繼解出

用matlab程式設計,編寫三次樣條函式(見附錄),對第一題求解:

>> format short g;

>> x1=[8.125,8.4,9.0,9.485,9.6,9.959,10.166,10.2]';

>> y1=[0.0774,0.099,0.280,0.60,0.708,1.200,1.800,2.177]';

>> u1=0.01087;un=100;

>> xx1=[x1(1):0.001:x1(end)]';

>> [yy1 b1 c1 d1]=spline3(x1,y1,xx1,1,u1,un);

>> fprintf('\t\tb1\t\tc1\t\td1\n');

b1 c1 d1

>> disp([b1 c1(1:end-1,1) d1]);

0.01087 0.14489 0.368

0.17405 0.4485 -0.393

0.2878 -0.25891 2.1153

1.5294 2.8188 -69.141

-0.56548 -21.035 73.614

12.794 58.247 -512.32

-28.949 -259.9 42279

>> plot(x1,y1,'bo',xx1,yy1,'r-');

>> grid on

畫出插值曲線的影象。

圖1 三次樣條曲線

2.逆時針旋轉座標軸45o 保持(1.)中結點和邊界條件的幾何關係不變,再次作三次樣條函式插值,畫出插值曲線的影象。

座標軸逆時針旋轉45°,相當於節點順時針旋轉45°。設為旋轉前的座標,為旋轉後的座標,則可以得到如下關係:

故旋轉後的節點座標為:

>> theta=-pi/4;

>> for i=1:length(x1)

x2(i)=cos(theta)*x1(i)-sin(theta)*y1(i);

y2(i)=sin(theta)*x1(i)+cos(theta)*y1(i);

end>>fprintf('\t\t\tx2\t\t\ty2\n');

>> disp([x2' y2']);

5.8 -5.6905

6.0097 -5.8697

6.562 -6.166

7.1312 -6.2826

7.2889 -6.2876

7.8906 -6.1935

8.4612 -5.9157

8.7519 -5.6731

端點處的一階導數為:

>> v1=(u1+tan(theta))/(1-u1*tan(theta));

>> vn=(un+tan(theta))/(1-un*tan(theta));

>> fprintf('\t\t\tv1\t\t\tvn\n');

v1vn

>> disp([v1 vn]);

-0.97849 0.9802

則旋轉後的三次樣條的係數及影象為:

>> xx2=[x2(1):0.001:x2(end)]';

>> [yy2 b2 c2 d2]=spline3(x2,y2,xx2,1,v1,vn);

>> fprintf('\t\t\tb2\t\t\tc2\t\t\td2\n');

b2c2d2

>> disp([b2 c2(1:end-1,1) d2]);

-0.97849 0.67221 -0.38277

-0.74704 0.43138 -0.090754

-0.35362 0.28102 -0.034909

-0.067629 0.22141 0.053338

0.0061747 0.24664 0.0046897

0.3081 0.2551 0.10233

0.6992 0.43028 0.12195

>> plot(x2,y2,'b+',xx2,yy2,'m-.');

>> grid on;

圖2 旋轉後的三次樣條曲線

3.比較(1.)、(2.)的結果,能得到什麼結論?

將(1.)中所得的三次樣條曲線整體順時針旋轉45°,並與二題(2.)中的三次樣條曲線畫在同一幅圖中比較,得

>> for i=1:length(xx1)

xx3(i)=cos(theta)*xx1(i)-sin(theta)*yy1(i);

yy3(i)=sin(theta)*xx1(i)+cos(theta)*yy1(i);

end>> plot(x2,y2,'bo',xx3,yy3,'r--',xx2,yy2,'m');

>> grid on;

>> legend('節點','旋轉前','旋轉後');

圖3旋轉前後三次樣條曲線幾何比較

比較圖中兩條曲線可知,曲線不重合,故三次樣條插值不具備幾何不變性。

附錄:三次樣條插值函式程式

function[yy,b,c,d]=spline3(x,y,xx,flag,vl,vr)

%三次樣條插值函式

%(x,y)為插值節點,xx為插值點

%flag表端點邊界條件型別;

%flag=0:自然樣條(端點二階導數為0);

%flag=1:第一類邊界條件(端點一階導數給定);

%flag=2:第二類邊界條件(端點二階導數給定);

%vl,vr表示左右端點處的在邊界條件值;

%樣條函式為:si(x)=yi+bi*(x-xi)+ci*(x-xi)^2+di*(x-xi)^3

%b,c,d分別為各子區間上的系數值

%yy表示插值點處的函式值

if length(x)~=length(y)

error('輸入資料應成對!');

endn=length(x);

a=zeros(n-1,1);

b=a;d=a;dx=a;dy=a;

a=zeros(n);b=zeros(n,1);

for i=1:n-1

a(i)=y(i);

dx(i)=x(i+1)-x(i);

dy(i)=y(i+1)-y(i);

endfor i=2:n-1

a(i,i-1)=dx(i-1);

a(i,i)=2*(dx(i-1)+dx(i));

a(i,i+1)=dx(i);

b(i,1)=3*(dy(i)/dx(i)-dy(i-1)/dx(i-1));

end%自然樣條端點條件(端點二階導數為0)

if flag==0

a(1,1)=1;

a(n,n)=1;

end%端點一階導數條件

if flag==1

a(1,1)=2*dx(1);a(1,2)=dx(1);

a(n,n-1)=dx(n-1);a(n,n)=2*dx(n-1);

b(1,1)=3*(dy(1)/dx(1)-vl);

b(n,1)=3*(vr-dy(n-1)/dx(n-1));

end%端點二階導數條件

if flag==2

a(1,1)=2;a(n,n)=2;

b(1,1)=vl;b(n,1)=vr;

endc=a\b;

for(i=1:n-1)

d(i)=(c(i+1)-c(i))/(3*dx(i));

b(i)=dy(i)/dx(i)-dx(i)*(2*c(i)+c(i+1))/3;

end[mm,nn]=size(xx);

yy=zeros(mm,nn);

for i=1:mm*nn

for ii=1:n-1

if xx(i)>=x(ii)&xx(i)j=ii;

break;

elseif xx(i)==x(n)

j=n-1;

endendyy(i)=a(j)+b(j)*(xx(i)-x(j))+c(j)*(xx(i)-x(j))^2+d(j)*(xx(i)-x(j))^3;

endend

計算方法第二次上機作業

gegebao 摘要 程式基於matlab,包括問題陳述 演算法與程式 結果與分析 討論四個部分。一 問題陳述 數學上已經證明了 成立,所以可以通過積分來計算的近似值。1 分別使用矩形 梯形和simpson復合求積公式計算的近似值。選擇不同的h,對於每種求積公式,試將誤差刻畫成h的函式,並比較各方面...

第二次作業

案例分析題 一 案情 四川省瀘州市某公司職工黃永彬和蔣倫芳1963年結婚,但是妻子蔣一直沒有生育,後來只得抱養了乙個兒子。由此原因給家庭籠罩上了一層陰影。1994年,黃永彬認識了乙個名叫張學英的女子,並且在與張認識後的第二年同居。黃的妻子蔣發現這一事實以後,進行勸告但是無效。1996年底,黃永彬和張...

圖論第二次作業

一 第四章 4.3 1 畫乙個有euler閉跡和hamilton圈的圖 2 畫乙個有euler閉跡但沒有hamilton圈的圖 3 畫乙個有hamilton圈但沒有euler閉跡的圖 4 畫乙個既沒有euler閉跡也沒有hamilton圈的圖 解 1 乙個有euler閉跡和hamilton圈的圖形如...