数值分析matlab程序
- 格式：doc
- 大小：785.50 KB
- 文档页数：46
数值分析（matlab程序）
曹德欣 曹璎珞
第一章 绪论
数值稳定性程序，计算P11 试验题一积分
function try_stable
global n
global a
a=input('a=');
N = 20;
I0 = log(1+a)-log(a);
I = zeros(N,1);
I(1) = -a*I0+1;
for k = 2:N
I(k) = - a*I(k-1)+1/k;
end
II = zeros(N,1);
if a>=N/(N+1)
II(N) = (1+2*a)/(2*a*(a+1)*(N+1));
else
II(N) =(1/((a+1)*(N+1))+1/N)/2;
end
for k = N:-1:2
II(k-1) = ( - II(k) +1/k) / a;
end
III = zeros(N,1);
for k = 1:N
n = k;
III(k) = quadl(@f,0,1);
end
clc
fprintf('\n 算法1结果 算法2结果 精确值')
for k = 1:N,
fprintf('\nI(%2.0f) %17.7f %17.7f %17.7f',k,I(k),II(k),III(k))
end
function y = f(x)
global n
global a
y = x.^n./(a+x);
return 第二章 非线性方程求解
下面均以方程y=x^4+2*x^2-x-3为例：
1、二分法
function y=erfen(a,b,esp)
format long
if nargin<3 esp=1.0e-4;
end
if fun(a)*fun(b)<0
n=1;
c=(a+b)/2;
while c>esp
if fun(a)*fun(c)<0
b=c;
c=(a+b)/2;
elseif fun(c)*fun(b)<0
a=c;
c=(a+b)/2;
else y=c; esp=10000;
end
n=n+1;
end
y=c;
elseif fun(a)==0
y=a;
elseif fun(b)==0
y=b;
else disp('these,nay not be a root in the intercal')
end
n
function y=fun(x)
y=x^4+2*x^2-x-3;
2、牛顿法
function y=newton(x0)
x1=x0-fun(x0)/dfun(x0);
n=1;
while (abs(x1-x0)>=1.0e-4) & (n<=100000000)
x0=x1;
x1=x0-fun(x0)/dfun(x0);
n=n+1;
end
y=x1
n
function y=fun(x) y=x^4+2*x^2-x-3;
3、割线法
function y=gexian(x0,x1)
x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %根据初始XO和X1求X2
n=1;
while (abs(x1-x0)>=1.0e-4) & (n<=100000000) %判断两个条件截止
x0=x1; %将x1赋给x0
x1=x2; %将x2赋给x1
x2=x1-fun(x1)*(x1-x0)/(fun(x1)-fun(x0)); %迭代运算
n=n+1;
end
y=x2
n
function y=fun(x)
y=x^4+2*x^2-x-3;
第四章
题目：推导外推样条公式：1232123211223322~~~~~22~nnnnnnnnddddMMMM，并编写程序与Matlab的Spline函数结果进行对比，最后调用追赶法解方程组。
解：设()sx的二阶导数值''()(0,1,,)iisxMin，因为()sx在1,iixx上是三次多项式，所以''()isx在1,iixx上是一次多项式，且可表示为
''11()iiiiiiixxxxsxMMhh （1）
对式（1）进行两次积分可得：
33111()()()()()66kkkkkkkkkkkMMsxxxxxPxxQxxhh （2）
由式（2）根据()kkysx和11()kkysx可得系数为： 6kkkkkymhPh，116kkkkkymhQh （3）
将式（3）代入式（2）可得()sx为：
3311111()()()()()()()6666kkkkkkkkkkkkkkkkkMMymhymhsxxxxxxxxxhhhh
分别对()kksx和1()kksx进行求导，并根据两者相等整理得：
112kkkkkkMMMd （4）
其中，11kkkkhhh，1kkkkhhh，11116()kkkkkkkkkyyyydhhhh
但是要解给定的方程组，还需要两个另外的条件，而外推样条插值的条件可通过下面推理得出：令式（1）中的x等于1x，可得：
32112122()MMMhhhhh （5）
1212122()()nnnnnnnnMMMhhhhh （6）
令式（4）中2k，然后将式（5）代入得：
1123222(2)(1)hhMMdhh （7）
令式（4）中1kn，然后将式（6）代入得：
1112122(2)(1)nnnnnnnhhMMdhh （8）
将方程组（4）和式（7）、式（8）联立展开即得题目所求。
按照推导出的外推样条插值公式编程可得如下M文件（spline_wt.m）
function spline_wt
X=[0 1 2 3 4 5];
Y=[0 0.5 2 1.5 3.5 1.9];
% 调用自编程序
S = spline_w(X,Y);
%调用matlab提供程序
S1=spline(X,Y);
fnplt(S,'r',2); % 作图
hold on
fnplt(S1,'b',1); hold on
plot(X,Y,'or'); % 画上节点
title('红线为自编程序曲线，蓝线为自带程序曲线')
% pp的第j行表示第j个三次多项式的4系数并写出分段多项式
pp = S.coefs
P1 = poly2str(pp(1,:),'(x-0)')
P2 = poly2str(pp(2,:),'(x-1)')
P3 = poly2str(pp(3,:),'(x-2)')
ppp=S1.coefs
PP1 = poly2str(ppp(1,:),'(x-0)')
PP2 = poly2str(ppp(2,:),'(x-1)')
PP3 = poly2str(ppp(3,:),'(x-2)')
function sp = spline_w(X,Y)n = length(X);
h = diff(X);
d = diff(Y)./h;
d1(2:n-1)=6*diff(d)./(h(1:n-2)+h(2:n-1));
mu(2:n-1)=h(1:n-2)./(h(1:n-2)+h(2:n-1));
mu(n-1)=1-h(n-2)/h(n-1);
la(2:n-1)=1-mu(2:n-1);
la(2)=1-h(1)/h(2);
% 计算三对角方程组
a = mu(3:n-1);
b = 2*ones(n-2,1);
b(1)=2+h(1)/h(2);
b(n-2)=2+h(n-2)/h(n-1);
c = la(2:n-2);
u(1:n-2) = d1(2:n-1);
% 调用追赶法解方程组
M(2:n-1) = tridiag(a,b,c,u); M(1)=(1+h(1)/h(2))*M(2)-M(3)*h(1)/h(2);
M(n)=(1+h(n-1)/h(n-2))*M(n-1)-M(n-2)*h(n-1)/h(n-2);
% 下面计算分段多项式的四个系数
S=zeros(n-1,4);
for k=0:n-2
S(k+1,1)=(M(k+2)-M(k+1))/(6*h(k+1));S(k+1,2)=M(k+1)/2;
S(k+1,3)=d(k+1)-h(k+1)*(2*M(k+1)+M(k+2))/6;
S(k+1,4)=Y(k+1);
end
sp = mkpp(X,S);%转换成 Matlab 格式
