數值計算第二次大作業
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圈的圖形如...