数值分析在MATLAB中的实现（M函数文件）

学习数值分析过程中编写的各章中涉及的方法的M函数
相互交流，如有错误请多指教
部分内容如下：（na——numerical analysis）
数据建模：
nalagr——拉格朗日插值
naspline——三阶样条插值（一阶导数边界条件）
nafit——多项式拟合

数值微积分：
natrapz——定步长梯形法
nagsint——定步长Gauss积分
naromberg——Romberg求积
naadapt——自适应步长Simpson法求积

常微分方程的数值解法：
naeuler——欧拉格式
naeuler2——两步欧拉格式，改进的欧拉格式
nark4——4阶龙格库塔格式
nark4v——变步长4阶龙格库塔格式

非线性方程及线性方程组的迭代解法：
nabisect——二分法（收敛较慢）
nanewton——Newton迭代法
nags——解普通方程组的Gauss-Seidel迭代
naspgs——解大型稀疏方程组的Gauss-Seidel迭代
nasor——分量形式的SOR迭代

线性方程组的直接解法：
nagauss——顺序Gauss消去法
nagauss2——选主元Gauss消去法
nalu——LU分解
nalupad——紧凑格式的LU分解

ma3biao1online

2010年04月26日 21:38:49

2楼

感谢各位捧场

[ 本帖最后由 ma3biao1online 于 2010-11-8 09:23 编辑 ]

ma3biao1online

2010年04月27日 13:00:37

3楼

function x=nasor(A,b,omega,x0,e,N)
%用途：用分量形式的SOR迭代解线性方程组Ax=b
%格式：x=nasor(A,b,omega,x0,e,N) A为系数矩阵，b为右端向量，x返回解向量，x0为初值向量（默认原点）
% e为精度（默认1e-4），设置迭代次数上限以防发散（默认500）,omega是松弛因子，一般取1-2之间的
% 数（默认1.5）
n=length(b);
if nargin<6,N=500;end %nargin为输入变量个数，此处指N缺省
if nargin<5,e=1e-4;end
if nargin<4,x0=zeros(n,1);end
if nargin<3,omega=1.5;end
x=x0;x0=x+2*e;
k=0;L=tril(A,-1);U=triu(A,1); %tril(A,-1)指严格下三角矩阵，triu(A,1)指严格上三角矩阵，
%tril(A,K)表示提取矩阵A第K条对角线上及其以下的各元素形成的矩阵,其中K表示第K条对角线，
%K=0，指主对角线，K>0,指在主对角线以上的第K条对角线，K<0则在主对角线以下
%triu(A,K)表示提取矩阵A第K条对角线上及其以上的各元素形成的矩阵
while norm(x0-x,inf)>e&k %1-范数，无穷大范数，p的缺省值为2 inf指无穷大
k=k+1,x0=x;
for i=1:n
x1(i)=(b(i)-L(i,1:i-1)*x(1:i-1,1)-U(i,i+1:n)*x0(i+1:n,1))/A(i,i);
x(i)=(1-omega)*x0(i)+omega*x1(i); %加权平均
end
disp(x') %输出矩阵x的转置
end
if k==N,warning('已达到迭代次数上限');end

[ 本帖最后由 ma3biao1online 于 2010-5-29 10:54 编辑 ]

ma3biao1online

2010年04月27日 13:05:54

4楼

function [ x,y] = nark4v( dyfun,xspan,y0,e,h )
%用途：变步长4阶经典Runge-Kutta格式解常微分方程y'=f(x,y),y(xo)=y0
%格式：[x,y]=nark4(defun,xspan,y0,h) dyfun为函数f(x,y),xspan为求解区间
% [x0,xN],y0为初值y（x0），h为初始步长（默认为xspan的1/10），x返回节点，
% y返回数值解，e为精度要求
%feval为函数求值命令，'为转置，nargin为输入变量个数,eps为浮点数识别精度2^（-52）
if nargin<5,h=(xspan(2)-xspan(1))/10;end
n=1;x(n)=xspan(1);y(n)=y0;
[y1,y2]=comput(dyfun,x(n),y(n),h);
while x(n) if abs(y2-y1)/10>e
while abs(y2-y1)/10>e
h=h/2;
[y1,y2]=comput(dyfun,x(n),y(n),h);
end
else
while abs(y2-y1)/10<=e
h=2*h;
[y1,y2]=comput(dyfun,x(n),y(n),h);
end
h=h/2;h=min(h,xspan(2)-x(n));
[y1,y2]=comput(dyfun,x(n),y(n),h);
end
n=n+1;x(n)=x(n-1)+h;y(n)=y2;
[y1,y2]=comput(dyfun,x(n),y(n),h);
end
x=x';y=y';
function [y1,y2]=comput(dyfun,x,y,h)
y1=rk4(dyfun,x,y,h);
y21=rk4(dyfun,x,y,h/2);
y2=rk4(dyfun,x+h/2,y21,h/2);
function y=rk4(dyfun,x,y,h)
k1=feval(dyfun,x,y);
k2=feval(dyfun,x+h/2,y+h/2*k1);
k3=feval(dyfun,x+h/2,y+h/2*k2);
k4=feval(dyfun,x+h,y+h*k3);
y=y+h*(k1+2*k2+2*k3+k4)/6;