Labfans是一个针对大学生、工程师和科研工作者的技术社区。 | 论坛首页 | 联系我们(Contact Us) |
![]() |
|
![]() |
#1 |
初级会员
注册日期: 2011-10-14
年龄: 36
帖子: 2
声望力: 0 ![]() |
![]()
x=[8.125 8.4 9.0 9.485 9.6 9.959 10.166 10.2];
y=[0.0774 0.099 0.280 0.60 0.708 1.200 1.800 2.177]; Y1=0.01087; Y2=100; n=size(x,2); for i=1:n-2 %形成系数矩阵 u(i)=(x(i+1)-x(i))/(x(i+2)-x(i)); end for i=1:n-2 la(i)=(x(i+2)-x(i+1))/(x(i+2)-x(i)); end A=2*eye(n-2,n-2)+diag(u(1:n-3),1)+diag(la(2:n-2),-1); for i=1:n-1 %形成非齐次项 f(i)=(y(i+1)-y(i))/(x(i+1)-x(i)); end for i=1:n-2 g(i)=3*(la(i)*f(i)+u(i)*f(i+1)); end g(1)=g(1)-la(1)*Y1; g(n-2)=g(n-2)-u(n-2)*Y2; disp(g); [m,L,U]=zhuigan_Y(A,g); m=[Y1;m';Y2]; syms x1 aerfa beita; syms s s1; for i=2:n aerfa(i)=((x1-x(i-1))^2/(x(i)-x(i-1))^2)*(1+2*(x1-x(i))/(x(i-1)-x(i))); end aerfa(1)=(x1-x(2))^2/(x(1)-x(2))^2*(1+2*(x1-x(1))/(x(2)-x(1))); for i=2:n beita(i)=((x1-x(i-1))/(x(i)-x(i-1)))^2*(x1-x(i)); end beita(1)=((x1-x(2))/(x(1)-x(2)))^2*(x1-x(1)); for i=1:n s(i)=y(i)*aerfa(i)+m(i)*beita(i); end for i=1:n-1 %样条插值在各个区间上的表达式 s1(i)=s(i)+s(i+1); end fplot(inline(s1(2)),[x(2),x(3)]); 着个区间上的图形有问题,第一类边界条件在左边界一阶倒数为Y1,右边界一阶导数为Y2 希望老鸟们帮助帮助小弟 |
![]() |
![]() |
![]() |
#2 |
初级会员
注册日期: 2011-10-14
年龄: 36
帖子: 2
声望力: 0 ![]() |
![]()
中间追赶法的程序
unction [x,L,U] = zhuigan_Y(A,g) % A:系数矩阵 % g:等号右端的向量 aerfa(1)=A(1,1); %产生追赶法的L,U矩阵 gama=diag(A,-1); c=diag(A,1); n=size(A,1); beita=zeros(1,n-1); for i=1:n-1 beita(i)=c(i)/aerfa(i); aerfa(i+1)=A(i+1,i+1)-gama(i)*beita(i); end L=diag(aerfa)+diag(gama,-1); U=eye(n,n)+diag(beita,1); y=zeros(1,n); y(1)=g(1)/A(1,1); for i=2:n y(i)=(g(i)-gama(i-1)*y(i-1))/(A(i,i)-gama(i-1)*beita(i-1)); end x(n)=y(n); for j=n-1:-1:1 x(j)=y(j)-beita(j)*x(j+1); end end |
![]() |
![]() |