数值分析中常用的matlab程序

1.

% 最小二乘法拟合数据点方法1：

% 左除右除：xA=B ==> x=B/A | Ax=B ==> x=B\A

% A 表示拟合函数的组合，如：多项式插值，A=[1,x,x.^2,...,x.^n],表示拟合函数为

% 多项式：s(x)=a0+a1*x+...+an*x^n;又如：A=[log(x),cos(x),exp(x)]

%则表示拟合函数为s(x)=a0*ln(x)+a1*cos(x)+a2*exp(x)

% 法方程为：A'*A*z=A'*y ==> A*z=y ==> z=A\y z=(a0,...,an)'

% Date: 2012-1-1

clear;

clc;

x=[0 0.25 0.5 0.75 1]';

y=[1 1.284 1.6487 2.1170 2.7183]';%索要拟合的数据点

x1=ones(size(x),1);

A=[x1 x x.^2];%拟合函数

Z=A\y %A中每列函数的参数

% 最小二乘法拟合数据方法2：采用polyfit函数

% Date:2012-1-1

clear;

clc;

x=[0 0.25 0.5 0.75 1]';

y=[1 1.284 1.6487 2.1170 2.7183]';%索要拟合的数据点

p=polyfit(x,y,2) %polyfit(x,y,n),(x,y)为数据点坐标，n为拟合多项式阶数，p 为

% 所求拟合多项式的幂次从高到低排列的系数。

2.

% 复合梯形公式计算积分值

% 输入：fun--积分函数；a，b--积分区间；n--区间等分数

% 输出：I--数值积分结果

% 调用格式(ex)：re=ftrapz(@fun1,0,1,10)

% 2012-1-1

function I=ftrapz(fun,a,b,n)

h=(b-a)/n;%区间等分

x=linspace(a,b,n+1);%将a到b的区间等分成(n+1)-1个区间，数据点有(n+!)个

y=feval(fun,x);

I=h*(0.5*y(1)+sum(y(2:n))+0.5*y(n+1));

%积分原函数

% Date：2012-1-1

function y=fun1(x)

y=exp(-x);

3.

% 复合simpson公式求积分

% 输入：fun--积分函数；a，b--积分区间；n--区间等分数% 输出：I--数值积分结果

% 调用格式(ex)：re=fsimpson(@fun1,0,1,10)

% 2012-1-1

function I=fsimpson(fun,a,b,n)

h=(b-a)/n;

x=linspace(a,b,2*n+1);

y=feval(fun,x);

I=(h/6)*(y(1)+2*sum(y(3:2:2*n-1))+4*sum(y(2:2:2*n))+y(2*n+1));

4.

% 两点GS-Legendre公式求积分

% 输入：fun--积分函数；a，b--积分区间；

% 输出：I--数值积分结果

% 调用格式(ex)：re=GSLege(@fun2,0,1)

% 2012-1-1

function I=GSLege(fun,a,b)

%将区间[a.b]通过变量替换x=(a+b)/2-(b-a)/2*t变到[-1,1]

%其中t取GS-Lege的高斯点

m1=feval(fun,(a+b)/2+(b-a)/2*(-1/sqrt(3)));

m2=feval(fun,(a+b)/2+(b-a)/2*(1/sqrt(3)));

I=(b-a)/2*(m1+m2);

% GS-Legendre公式的积分函数

% Date:2012-1-1

function y=fun2(x)

y=sin(x)/x;

5.

% 追赶法求解线性方程组Ax=b，其中A是三对角方阵

%function x=tridiagsolver(A,b)

clear;

clc;

A=[2 -1 0 0;-1 3 -2 0;0 -2 4 -3;0 0 -3 5];%三对角矩阵，线性方程组系数矩阵b=[6,1,-2,1]';%

[n,n]=size(A);

for i=1:n

if(i<2)

l(i)=A(i,i);

y(i)=b(i)/l(i);

u(i)=A(i,i+1)/l(i);

elseif i

l(i)=A(i,i)-A(i,i-1)*u(i-1);

y(i)=(b(i)-y(i-1)*A(i,i-1))/l(i);

u(i)=A(i,i+1)/l(i);

else

l(i)=A(i,i)-A(i,i-1)*u(i-1);

y(i)=(b(i)-y(i-1)*A(i,i-1))/l(i);

end

end

x(n)=y(n);

for j=n-1:-1:1

x(j)=y(j)-u(j)*x(j+1);

end

x

6.

% SOR迭代求解非线性方程组Ax=b

% 输入：A--系数矩阵；b--；omega--松弛因子（0~2);tol--精度

% 输出：x--方程的解向量;iter--迭代次数;

% 调用格式(ex)：[x,iter=sor(A,b,1.1,1e-4)

% 2012-1-1

%function [x,iter]=sor(A,b,omega,tol)

%{%}

clear;clc;

A=[2 -1 0;-1 3 -1;0 -1 2];

b=[1 8 -5]';

omega=1.1;

tol=1e-4;

D=diag(diag(A));%diag(A)返回的是A的对角元组成的列向量；diag(b)返回的是以列向量b为对角元的方阵；

L=D-tril(A);%tril(A)返回的是矩阵A的下三角矩阵；

U=D-triu(A);%triu(A)返回的是矩阵A的上三角矩阵；

x=zeros(size(b));%给定迭代初始值零向量

iter=1;

while iter<500

x=(D-omega*L)\(omega*b+(1-omega)*D*x+omega*U*x);%迭代格式

error=norm(b-A*x)/norm(b);%收敛条件

if error

break;

end

iter=iter+1;

end

if iter>= 500

fprintf('root not found!');

end

7.

% 牛顿法求解非线性方程的根

% 输入：fun--非线性函数；dfun--非线性函数导数；x0--初始值；tol--精度；

% 输出：x--非线性方程数值根

% 调用格式(ex)：x=newton(@fun3,@dfun3,6,1e-3)

% 2012-1-1

function x=newton(fun,dfun,x0,tol)

iter=1;

if abs(feval(fun,x0))

x=x0;

return;

end

while iter<500

if iter==1

x=x0-feval(fun,x0)/feval(dfun,x0);

if abs(feval(fun,x))

break;

end

else

x=x-feval(fun,x)/feval(dfun,x);

if abs(feval(fun,x))

break;

end

end

iter=iter+1;

end

if iter==500

fprintf('not successful!');

x=NaN;

end

% newton的函数文件

% Date：2012-1-1

function y=fun3(x)

y=3.*x.^3-8.*x.^2-8.*x-11;

% newton的导函数文件

% Date：2012-1-1

function y=dfun3(x)

y=9.*x.^2-16.*x-8;

8.

% 两点割线法求解非线性方程的根

% 输入：fun--非线性函数；a,b--两个初始值；tol--精度；

% 输出：x--非线性方程数值根

% 调用格式(ex)：x=gexian(@fun3,3,6,1e-3)

% 2012-1-1

function x=gexian(fun,a,b,tol)

iter=1;

xk1=a;

xk2=b;

while iter<500

xk3=xk2-feval(fun,xk2)*(xk2-xk1)/(feval(fun,xk2)-feval(fun,xk1));

if abs(feval(fun,xk3))

break;

else

xk1=xk2;

xk2=xk3;

end

iter=iter+1;

end

if iter==500

fprintf('not successful!');

x=NaN;

else

x=xk3;

end

9.

% 乘幂法求矩阵的按模最大特征值及其特征向量

% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；% 输出：x--特征向量；lam--按模最大特征值

% 调用格式(ex)：[lam,x]=power(A,v0,1e-2)

% 2012-1-1

function [lam,x]=matrixpower(A,v0,tol)

v1=A*v0;

v2=A*v1;

sum=0;

p=0;

for i=1:size(v1)

if v1(i)~=0

sum=sum+v2(i)/v1(i);

p=p+1;

end

end

lam0=sum/p;

iter=2;

vk1=v2;

while iter<500

vk2=A*vk1;

sum=0;

p=0;

for i=1:size(vk1)

if vk1(i)~=0

sum=sum+vk2(i)/vk1(i);

p=p+1;

end

end

lam=sum/p;

if abs(lam-lam0)

break;

else

lam0=lam;

vk1=vk2;

end

end

if iter==500

fprintf('not successful!');

lam=NaN;

else

x=vk2;

end

9_2

% 改进乘幂法求矩阵的按模最大特征值及其特征向量

% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；

% 输出：x--特征向量；lam--按模最大特征值

% 调用格式(ex)：[lam,x]=pMatrixPower(A,v0,1e-2)

% Date:2012-1-2

function [lam,x]=pMatrixPower(A,v0,tol)

[ty,ti]=max(abs(v0));%返回v0中元素绝对值最大的元素值与下标,ti为下标lam0=v0(ti);

u0=v0/lam0;

iter=1;

while iter<500

v1=A*u0;

[tv,ti]=max(abs(v1));

lam1=v1(ti);

u0=v1/lam1;

if abs(lam0-lam1)

break;

else

lam0=lam1;

iter=iter+1;

end

end

if iter>=500

fprintf('not successful!');

lam=NaN;

else

lam=lam1;

x=u0;

end

10.

% 反幂法求矩阵的按模最小特征值及其特征向量

% 输入：A--要求的矩阵；v0--初始非零向量；tol--精度；

% 输出：x--特征向量；lam--按模最大特征值

% 调用格式(ex)：[lam,x]=invMatrixPower(A,v0,1e-2)

% Date:2012-1-2

function [lam,x]=invMatrixPower(A,v0,tol)

[ty,ti]=max(abs(v0));

lam0=v0(ti);

u0=v0/lam0;

iter=1;

while iter<500

v1=A\u0;

[tv,ti]=max(abs(v1));

lam1=v1(ti);

u0=v1/lam1;

if abs(1/lam0-1/lam1)

break;

else

lam0=lam1;

iter=iter+1;

end

end

if iter>=500

fprintf('not successful!');

lam=NaN;

else

lam=1/lam1;

x=u0;

end

11.

% 欧拉方法求解一阶常微分方程初值问题

% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;

% 输出：y--常微分方程在区间[a,b]上各点的数值解；

% 调用格式(ex)：y=odeEuler(@fun4,0,1,1,10)

% Date:2012-1-2

function y=odeEuler(fun,a,b,y0,n)

h=(b-a)/n;

y(1)=y0;

x=a:h:b;

for i=1:n

y(i+1)=y(i)+h*feval(fun,x(i),y(i));

end

% 常微分方程fun4=0

% Date：2012-1-2

function re=fun4(x,y)

re=y-2*x/y;

12.

% 改进欧拉公式（预估--校正）方法求解一阶常微分方程初值问题

% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;

% 输出：y--常微分方程在区间[a,b]上各点的数值解；

% 调用格式(ex)：y=pOdeEuler(@fun4,0,1,1,10)

% Date:2012-1-2

function y=pOdeEuler(fun,a,b,y0,n)

h=(b-a)/n;

y(1)=y0;

x=a:h:b; for i=1:n

yp=y(i)+h*feval(fun,x(i),y(i));

yc=y(i)+h*feval(fun,x(i+1),yp);

y(i+1)=0.5*(yp+yc);

end

13.

% 梯形方法求解一阶常微分方程初值问题

% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;

% 输出：y--常微分方程在区间[a,b]上各点的数值解；

% 调用格式(ex)：y=trapezium(@fun4,0,1,1,10)

% Date:2012-1-2

function y=trapezium(fun,a,b,y0,n)

h=(b-a)/n;

y(1)=y0;

x=a:h:b;

tol=1e-6;

for i=1:n

%用不动点迭代的方法求解非线性方程：

%y(i+1)=y(i)+h*feval(fun,x(i),y(i))/2+h*feval(fun,x(i),y(i+1))/2;

iter=1;

yy0=1+(i-1)*h;%迭代初始值

while iter<500

yy1=y(i)+h*feval(fun,x(i),y(i))/2+h*feval(fun,x(i)+h,yy0)/2;

if abs(yy1-yy0)

break;

else

yy0=yy1;

iter=iter+1;

end

end

if iter>=500

printf('not successful!');

y=NaN;

return;

else

y(i+1)=yy1;

end

end

14.

% 标准四阶四段龙格-库塔方法求解一阶常微分方程初值问题

% 输入：fun--一阶常微分函数；a,b--求解区间；y0--函数在a点值y(a)；n--所分区间数;

% 输出：y--常微分方程在区间[a,b]上各点的数值解；

% 调用格式(ex)：y=longekuta(@fun4,0,1,1,10)

% Date:2012-1-2

function y=longekuta(fun,a,b,y0,n)

h=(b-a)/n;

y(1)=y0;

for k=2:n+1

x=a+(k-2)*h;

k1=h*feval(fun,x,y(k-1));

k2=h*feval(fun,x+h/2,y(k-1)+k1/2);

k3=h*feval(fun,x+h/2,y(k-1)+k2/2);

k4=h*feval(fun,x+h,y(k-1)+k3);

y(k)=y(k-1)+(k1+2*k2+2*k3+k4)/6;

end

15.插值

% 拉格朗日插值

% 输入：x,y--插值数据点（x，y均为行向量）；xh--要插值的点; % 输出：yh--插值结果；

% 调用格式(ex)：yh=lagrange(x,y,xh)

% Date:2012-1-2

function yh=lagrange(x,y,xh)

n=length(x);

m=length(xh);

yh=zeros(1,m);

c1=ones(n-1,1);

c2=ones(1,m);

for i=1:n

xp=x([1:i-1 i+1:n]);

yh=yh+y(i)*prod((c1*xh-xp'*c2)./(x(i)-xp'*c2));

%prod对输入的一个向量返回其所有分量的乘积

end

%拉格朗日调用

clear;

clc;

x=[11 12];

y=[2.3979 2.4849];

xh=11.75;

yh=lagrange(x,y,xh)

15_2

% 牛顿插值

% 输入：x,y--插值数据点（x，y均为行向量）；xh--要插值的点; % 输出：yh--插值结果；

% 调用格式(ex)：yh=newtonPol(x,y,xh)

% Date:2012-1-2

function yh=newtonPol(x,y,xh)

n=length(x);

p(:,1)=x;

p(:,2)=y;

for j=3:n+1

p(1:n+2-j,j)=diff(p(1:n+3-j,j-1))./(x(j-1:n)-x(1:n+2-j))';%求差商表end

q=p(1,2:n+1)';%求牛顿法的系数--取第一行

yh=0;

m=1;

yh=q(1);

for i=2:n

m=q(i);

for j=2:i

m=m*(xh-x(j-1));%求牛顿法中各多项式值(xh-x0)…(xh-x n-1) end

yh=yh+m;%求和

end

%牛顿插值调用

clear;

clc;

x=[11 12 13];

y=[2.3979 2.4849 2.5649];

xh=11.75;

yh=newtonPol(x,y,xh)

Matlab中常用的数值计算方法

Matlab中常用的数值计算方法 数值计算是现代科学和工程领域中的一个重要问题。Matlab是一种用于数值计算和科学计算的高级编程语言和环境，具有强大的数值计算功能。本文将介绍Matlab中常用的数值计算方法，包括数值积分、数值解微分方程、非线性方程求解和线性方程组求解等。 一、数值积分 数值积分是通过数值方法来近似计算函数的定积分。在Matlab中，常用的数值积分函数是'quad'和'quadl'。'quad'函数可以用于计算定积分，而'quadl'函数可以用于计算无穷积分。 下面是一个使用'quad'函数计算定积分的例子。假设我们想计算函数f(x) = x^2在区间[0, 1]上的定积分。我们可以使用如下的Matlab代码： ``` f = @(x) x^2; integral = quad(f, 0, 1); disp(integral); ``` 运行这段代码后，我们可以得到定积分的近似值，即1/3。 二、数值解微分方程 微分方程是描述自然界各种变化规律的数学方程。在科学研究和工程应用中，常常需要求解微分方程的数值解。在Matlab中，可以使用'ode45'函数来求解常微分方程的数值解。'ode45'函数是采用基于Runge-Kutta方法的一种数值解法。

下面是一个使用'ode45'函数求解常微分方程的例子。假设我们想求解一阶常微分方程dy/dx = 2*x，初始条件为y(0) = 1。我们可以使用如下的Matlab代码：``` fun = @(x, y) 2*x; [x, y] = ode45(fun, [0, 1], 1); plot(x, y); ``` 运行这段代码后，我们可以得到微分方程的数值解，并绘制其图像。 三、非线性方程求解 非线性方程是指方程中包含非线性项的方程。在很多实际问题中，我们需要求解非线性方程的根。在Matlab中，可以使用'fsolve'函数来求解非线性方程的根。 下面是一个使用'fsolve'函数求解非线性方程的例子。假设我们想求解方程x^2 - 2 = 0的根。我们可以使用如下的Matlab代码： ``` fun = @(x) x^2 - 2; x = fsolve(fun, 1); disp(x); ``` 运行这段代码后，我们可以得到方程的近似根，即约等于1.4142。 四、线性方程组求解

Matlab100个实例程序

程序代码： (代码标记[code]...[/code]) 1-32是：图形应用篇 33-66是：界面设计篇 67-84是：图形处理篇 85-100是：数值分析篇 实例1：三角函数曲线（1） function shili01 h0=figure('toolbar','none',... 'position',[198 56 350 300],... 'name','实例01'); h1=axes('parent',h0,... 'visible','off'); x=-pi::pi; y=sin(x); plot(x,y); xlabel('自变量X'); ylabel('函数值Y'); title('SIN( )函数曲线'); grid on 实例2：三角函数曲线（2） function shili02 h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例02'); x=-pi::pi; y=sin(x)+cos(x); plot(x,y,'-*r','linewidth',1); grid on xlabel('自变量X'); ylabel('函数值Y'); title('三角函数');

实例3：图形的叠加 function shili03 h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例03'); x=-pi::pi; y1=sin(x); y2=cos(x); plot(x,y1,... '-*r',... x,y2,... '--og'); grid on xlabel('自变量X'); ylabel('函数值Y'); title('三角函数'); 实例4：双y轴图形的绘制 function shili04 h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例04'); x=0:900;a=1000;b=; y1=2*x; y2=cos(b*x); [haxes,hline1,hline2]=plotyy(x,y1,x,y2,'semilogy','plot'); axes(haxes(1)) ylabel('semilog plot'); axes(haxes(2)) ylabel('linear plot'); 实例5：单个轴窗口显示多个图形 function shili05 h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例05'); t=0:pi/10:2*pi;

数值分析Matlab作业

1.对方程f(x)=0，先用二分法得到一个较好的初始解，再用牛顿迭代法求一个精度较高的解。 一、Matlab编程 1.先用二分法得到较好初始解 syms x; fun=input('(输入函数形式)fx='); a=input('（输入二分法下限）a='); b=input('（输入二分法上限）b='); d=input('输入误差限d=')%二分法求根 %f=inline(x^2-4*x+4); %修改需要求解的inline函数的函数体 f=inline(fun);%修改需要求解的inline函数的函数体 e=b-a; k=0 ; while e>d c=(a+b)/2; if f(a)*f(c)<0 b=c; elseif f(a)*f(c)>0 a=c; else a=c;b=c end e=e/2; k=k+1; end x=(a+b)/2; x%x为答案 k%k为次数 2.初始解作为初值，利用牛顿迭代法求精度较高的解 function [x,k]=Newtondd(f,x0,e) %%牛顿迭代法，求f(x)=0在某个范围内的根。 %%f为f(x)，x0为迭代初值（选用上述初始解），e为迭代精度，k为迭代次数 x_a=x0; x_b=x_a-subs(f,x_a)/subs(diff(f),x_a); k=1; while abs(x_a-x_b)>e, k=k+1; x_a=x_b; x_b=x_a-subs(f,x_a)/subs(diff(f),x_a); end

x=x_b; （x_b即为要求的精度更高的解）

二、实例验证 用二分法和牛顿迭代法分别计算方程sinx-2 x 2 =0的实根，精度为0.0001。 1.二分法求解较好的初始解 首先画出函数的关系图，可知x1=0,x2在1-1.5之间，即只需求解x2。 用二分法如上程序计算方程的根，初始区间为[1，1.5]。可得x2=1.4044，迭代次数为13次。 2.以上述初始解作为牛顿迭代法的初值，求解更高精度的解 取初值x0=1.4044，计算高精度解，可得x=1.4058，迭代次数为5次。 2.对非奇异矩阵A ，做A 的Doolittle 三角分解。 （最好用紧凑形式，如有困难 也可作通常的三角分解） 一、Matlab 编程 clear all clc A=input('请输入一个方阵 ');%输入一个n 阶方阵 [n,n]=size(A); L=zeros(n,n); U=zeros(n,n); for i=1:n %将L 的主对角线元素赋值1 L(i,i)=1; end for j=1:n %求矩阵U 的第一行元素 U(1,j)=A(1,j); end for k=2:n %求矩阵L 的第一列元素 L(k,1)=A(k,1)/U(1,1); end for i=2:n %求L 、U 矩阵元素 for j=i:n s=0; for t=1:i-1 s=s+L(i,t)*U(t,j); end U(i,j)=A(i,j)-s; end for k=i+1:n r=0; for t=1:i-1 r=r+L(k,t)*U(t,i); end L(k,i)=(A(k,i)-r)/U(i,i); end end

数值计算与金融仿真matlab程序总结(双面打印)

变量合法命名：只含字母、数字、下划线，并以字母开头 who 显示定义变量 whos 显示定义变量及其信息 lookfor “?”查询与关键词有关的所有函数 Integer 整数 real 实数 complex 复数 inf 无限大 NaN 非数字 format long 显示多位小数 short 4位 bank 2位 round()四舍五入 fix()去小数部分 prod() 连乘 asin() transpose()或矩阵 矩阵转置 eye(4)4行4列单位矩阵 q13(2,4) 矩阵第二行第四列数据 q13(2:3,2:3) 矩阵q13 2到3行，2到3列的数据 .*对应元素相乘 同理./ .^ a ： b ： c 以a 开始，以b 为间隔，最大数小于等于c 的数列 sum(A)对列求和(得行) sum(A ，2)对行 sum(sum(A))对列求和再对行求 prod(A)对列连乘 min(A)输出每一列的最小值 [m,j]=max(B)，m=最大值取值，j=最大值位置 判断：a==b a 等于b a~=b a 不等于b 正确输出1 . plot(x,y, '(颜色)(形状)') linspace(a,b,n)把起点为a 终点为b 的直线等分为n 份 plot(x,y,'.',x,cos(3*pi*x),'g*')两条线一图 legend('Sin curve', 'Cos curve')加图标 grid 加格子 hold on 保留原曲线，可用于画多线在一图 subplot(121), plot(x,y) 将作图区域分为1行2列，作第1个区域的图 N = 100; h = 1/N; x = 0:h:1; y = cos(3*pi*x); plot(x,y) x = linspace (0,1,101); y = sin(3*pi*x); plot(x,y) rand(3)3行方阵随机矩阵 randn r6 =random('Poisson',6,1,6) 文档读写IBM = xlsread('IBM.xls');IBM = textread('IBM.tx t');IBM = dlmread('IBM.txt','',1,1); IBM = importdata('IBM.xls');m 文件which mfile(文件名) 查看路径edit mfile 编辑脚本mfile 与函数mfile 的区别，前者是全局变量，后者是局部变量 function [var1 var2 …] = functionname(arg1,arg2, …) 多个图形x = -1:.05:1; for n = 1:8 subplot(4,2,n); plot(x,sin(n*pi*x)); end 求和s = 0 for i = 1:100 s = s + i; end sum 能以1^2 + 2^2 + … + n^2表示并小于100的? S = 1; n = 1; while S+ (n+1)^2 < 100 n = n+1; S = S + n^2; end [n, S] If...Then Loop S = 1; n = 1; for i=1:100；if S+(n+1)^2 < 100 S = S+ (n+1)^2; n = n+1; end end [n, S] Chapter 2 现金流分析 pvvar 求变动现金流量的现值pvfix fvvar fvfix PV = pvvar(CashFlow, Rate, CFDates) > CashFlow = [100 300 450] CFDates = ['01/12/1987'; '02/14/1988'; '03/03/1988'] PV = pvvar(CashFlow, 0.09, CFDates) 日期列数据 >> pvvar([-15000 3000 4500 5000 6800], 0.08) ans = 603.1667 没有现金流日期时，默认相隔一年 Irr 固定周期的内部收益率X irr 变动周期的内部收益率 Return = irr(CashFlow) Return = xirr(CashFlow, Dates) >CF = [-10000 3000 4500 5000]; Return = irr(CF) >>Return =0.1105 >Dates=['01/12/00'; '02/14/01'; '09/03/01'; '12/31/02']; Return = xirr(CF, Dates) >>Return =0.1177 cfdur and cfconv 久期和凸性 >CashFlow=[5 5 5 5 5 105];> [Dur, ModDur] = cfdur(CashFlow, 0.05)收益 Dur =5.3295ModDur =5.0757 > Conv = cfconv(CF, 0.05) Conv =95.5410 债券价格和收益率 [P, I 应付利息] = bndprice(Yield, CouponRate 票面利率, …Settle ?交收日, …Maturity ?到期) p=price+accruedint Yield = bndyield(Price, CouponRate, Settle, Maturity) Duration = bnddurp(Price, CouponRate, Settle, Maturity) Duration = bnddury(Yield, CouponRate, Settle, Maturity) Convex ity = bndconvp(Price, CouponRate, Settle, Maturity) p-y Price-yield curve >> yields = 0.01: 0.01: 0.20;>> [P I] = bndprice(yields, 0.1, '08/10/07', '12/31/20');>> plot (yields, P+I);>> grid on;>> xlabel('Yield');>> ylabel('Price') ;>> Title('Price-Yield Curve'); 利率免疫 duration=holding period >> settle = '28-Aug-2007';>> maturity = ['15-Jun-2012' ; '31-Oct-2017' ; '01-Mar-2027'];>> couponRate = [0.07 ; 0.06 ; 0.08];>> yield = [0.06 ; 0.07 ; 0.075];>> duration = bnddury(yield, couponRate, settle, maturity); >> convex ity = bndconvy(yield, couponRate, settle, maturity);% COMPUT E PORTFOLIO WEIGHTS >> A = [duration'; convexity'; 1 1 1];>> b = [ 10; 160; 1];>> weights = A\b 有效前沿 [PortRisk, PortReturn, PortWts] = frontcon(ExpReturn 历史收益率, ExpCovariance 历史协方差, Num Ports 返回结果的个数, PortReturn 目标收益, AssetBounds 资产界限, Groups, GroupBounds 组合界限) r = [0.2 0.1];两种资产历史收益率> s = [0.2 -0.1; -0.1, 0.4];历史协方差 >> [Risk, Return, Wts] = frontcon(r, s, 5);>> [Risk, Return, Wts] ans =0.2958 0.1625 0.6250 0.3750 0.3075 0.1719 0.7187 0.2813 0.3400 0.1812 0.8125 0.1875 0.3883 0.1906 0.9062 0.0938 0.4472 0.2000 1.0000 0 最后两列为权重 含约束条件的有效前沿 [PortRisk, PortReturn, PortWts] = portopt (Return, Cov, [], PortReturn, ConSet) 约束条件ConSet = portcons ('Default', Num, 'AssetLims', Min, Max,'GroupLims', Group, GroupMin, GroupMax); >> r = [0.1 0.2 0.15]; >> s = [0.0100,-0.0061,0.0042; -0.0061, 0.0400,-0.0252; 0.0042,-0.0252,0.0225];>> ObjReturn = 0.17; >> [Risk, Return, Wts] = portopt(r, s, [], ObjReturn) Return =0.1700 Wts =0.0278 0.4278 0.5445 >> Return = [0.1 0.2 0.15]; >> Cov = 0.01*[0.5,-1,0.4; -1,4, -0.2; 0.4,-0.2,2.3]; >> PortReturn = [0.15 0.16];目标收益率在0.15~0.16之间>> Num = 3; >> Min = [0.20 NaN NaN];投资组合中每种资产的下限>> Max = [0.5 0.5 0.5]; ……上限>> Group = [1 1 0; 0 1 1];分为两组，第一组由第1、2种资产构成；第二组……>> GroupMin = [0.2, 0.2]; 第一组的投资比例下限为0.2；二…… >> GroupMax = [0.8, 0.8]; 第一组的投资比例上限为0.8；二…… >> ConSet = portcons('Default', Num, 'AssetLims', Min, Max,'GroupLims', Group, GroupMin, GroupMax); >> [PortRisk, PortReturn, PortWts] = portopt(Return, Cov, [], PortReturn, ConSet) PortRisk = 0.0724; 0.0910 PortReturn = 0.1500; 0.1600 PortWts = 0.4000 0.4000 0.2000 0.2606 0.4606 0.2789 可借贷无风险资产时的有效前沿 [RiskyRisk, RiskyReturn, RiskyWts, RiskyFraction 所有风险资产占总资产的比例, OverallRisk, OverallReturn] = portalloc(PortRisk, PortReturn, PortWts, RisklessRate, BorrowRate, RiskAversion 风险厌恶程度) >> Ret = [0.1 0.2 0.15]; >> Cov = [ 0.005 -0.010 0.004; -0.010 0.040 -0.002; 0.004 -0.002 0.023]; >> [Risk, Return, Wts] = frontcon(Ret, Cov, 20); >> RiskfreeR = 0.08; BorrowR = 0.12; >> RiskA version = 3; >> portalloc (Risk, Return, Wts, RiskfreeR, BorrowR, RiskA version); (return the graphic) >>[RiskyRisk, RiskyReturn, RiskyWts, RiskyFraction, OverallRisk, OverallReturn] = portalloc (Risk, Return, Wts, RiskfreeR, BorrowR, RiskA version) 最优化问题 [x, fval] = fmincon(fun,目标函数, x0,从x0开始试值, A, b, Aeq, beq, lb, ub, nonlcon 非线性约束,options) options = optimset('LargeScale','off') 注意：x1要加括号x(1)；目标函数要建立m 文件；非线性约束要建立m 文件 程序 1: function [f] = Myobj(x) f=exp(x(1))*(4*x(1)^2+2*x(2)^2+4*x(1)*x(2)+x(2)+1); 2: 写约束条件function [c, ceq] = Mycon(x)% Nonlinear inequality constraints c=[1.5+x(1)*x(2)-x(1)-x(2); -x(1)*x(2)-10];% Nonlinear equality constraints ceq = []; Step 3:调用最优化路径：A = [1 1; 2 -1]; b = [20; 10];Aeq = [];beq = []; lb = [-10 -10];ub = [10, 10];x0 = [-1,1]; % Make a starting guess at the solution options = optimset('LargeScale','off'); [x, fval] = fmincon(@Myobj, x0, A, b, Aeq, beq, lb, ub, @Mycon, options) B-S 模型 [Call, Put] = blsprice(Price, Strike, Rate, T ime, V olatility, Yield) [CallDelta, PutDelta] = blsdelta(Price, Strike, Rate, T ime, V olatility, Yield) [CallEl, PutEl] = blslambda,[CallRho, PutRho]= blsrho, V olatility = blsimpv(Price, Strike, Rate, Time, V alue, Limit, Yield, T olerance, Class) Chapter 3 二叉树模型 European and American Call European Put American Put 欧式二叉树 function [] = BinomialEuro() T ,n,K,r,sigma,S0;deltat=T/n; u=exp(sigma*sqrt(deltat)); d=exp(-sigma*sqrt(deltat)); p=(exp(r*deltat)-d)/(u-d); q=1-p; for i=1:n+1 s(i)=u^(n+1-i)*d^(i-1)*S0; 股价 w(i)=nchoosek(n,i-1)*p^(n+1-i)*q^(i-1); x(i)=max(s(i)-K,0); y(i)=max(-s(i)+K,0); cv(i)=w(i)*x(i); pv(i)=w(i)*y(i); end C=sum(cv(:))*exp(-r*T) P=sum(pv(:))*ex p(-r*T) 美式二叉树 function [c]=BinomialAmer() deltat=T /n; u=exp(sigma*sqrt(deltat)); d=exp(-sigma*sqrt(deltat)); p=(exp(r*deltat)-d)/(u-d); q=1-p; s(1,1)=S0; for t=2:n+1 (r 扣除红利) for i=1:t s(t,i)=u^(t-i)*d^(i-1)*S0; end; end for i=1:n+1; x(n+1,i)=max(K-s(n+1,i),0); y(n+1,i)=max(s(n+1,i)-K,0); end for t=n:(-1):1 for i=1:t x(t,i)=(p*x(t+1,i)+q*x(t+1,i+1)) /exp(r*deltat); x(t,i)=max(x(t,i),K-s(t,i)); y(t,i)=(p*y(t+1,i)+q*y(t+1,i+1))/exp(r*deltat); y(t,i)=max(y(t,i),-K+s(t,i)); end; end c=y(1,1); p=x(1,1) Chapter 4 蒙特卡罗模拟 计算面积： y = cos x [-pi/2, pi/2] function [area] = MC1(n) X=rand(n,1)*pi-pi/2；Y=rand(n,1) ；Z=cos(X); counter=0; for i = 1 : n i f Y(i) <= Z(i) counter=counter+1; end end area=pi*counter/n; 估计pi function [Mypi] = Mypi() n = 10000；count = 0; for i = 1:n x = rand; y = rand; i f x^2+y^2 <= 1 count = count + 1; end end Mypi = 4*count/n; 计算积分 count = 0; for i = 1:n ；x = 2*rand+1; count = count + cos(x)*sin(x); end Myinte = (b-a)*count/n; 欧式看涨期权定价 Chapter 5 随机数生成器 线性同余产生器 x(i+1)=ax(i)mod m; u(i+1)=x(i+1)/m f unction[]=line(a,m,n) x(1)=2 %余数mod(被除数，除数) f or i=1:n x(i+1)=mod(a*x(i),m); u(i+1)=x(i+1)/m; end x u 混合线性同余产生器 x(i+1)=(ax(i)+b)mod m; u(i+1)=x(i+1)/m f unction [u] = lcg(seed, n) a = 1229;b = 1;m = 32768; x = zeros(1,n); x(1)=mod(a*seed+b, m); f or i = 2:n x(i)= mod(a*x(i-1)+b, m); end u = x./m; % make the random unmber between 0 and 1 plot(u(1:n-1), u(2:n), '.'); % plot the pair of point (u(n), u(n+1)) title('Autocorrelaton of LCG'); 计算循环的周期 f unction [p]=lcg(a,b,m,x0,n) x(1)=mod((a*x0+b),m) 1 22 1122121212121212: m in (4241) ..: 1.50 100 20 210 -10,10 x x o b j e x x x x x s t x x x x x x x x x x x x ++++--+≤--≤+≤-≤≤≤(),,r q t t t e d p u e d e u d σ σ -??-?-? = ==-0,0,,0 max{,0)T rT T i T i i f e p S K -==-∑ 0,0,,0 max{,0)T rT T i T i i f e p K S -==-∑ ,01,1,1max{max{,0}, ()}t i i r t t i t i t i f K S u d e pf qf --?+++=-+1 ()[()]()[()]()b b N i i a a b a f x dx E f x dx b a E f x f x N =-≈=-≈ 1()[()]()[()]()b b N i i a a b a f x dx E f x dx b a E f x f x N =-≈=-≈∑??()() ()221 1,, ()(0) ()? i i r T T Z i rT i i n n for i n generate Z set S T S e set c e S T K set c c c n σ-++ -===-=+

Matlab100个实例程序

程序代码：(代码标记[code]...[/code] ) 1-32是：图形应用篇 33-66是：界面设计篇 67-84是：图形处理篇 85-100是：数值分析篇 实例1：三角函数曲线（1） function shili01 h0=figure('toolbar','none',... 'position',[198 56 350 300],... 'name','实例01'); h1=axes('parent',h0,... 'visible','off'); x=-pi:0.05:pi; y=sin(x); plot(x,y); xlabel('自变量X'); ylabel('函数值Y'); title('SIN( )函数曲线'); grid on 实例2：三角函数曲线（2） function shili02 h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例02'); x=-pi:0.05:pi; y=sin(x)+cos(x); plot(x,y,'-*r','linewidth',1); grid on xlabel('自变量X'); ylabel('函数值Y'); title('三角函数');

实例3：图形的叠加 function shili03 h0=figure('toolbar','none',... 'position',[200 150 450 350],... 'name','实例03'); x=-pi:0.05:pi; y1=sin(x); y2=cos(x); plot(x,y1,... '-*r',... x,y2,... '--og'); grid on xlabel('自变量X'); ylabel('函数值Y'); title('三角函数'); 实例4：双y轴图形的绘制 function shili04 h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例04'); x=0:900;a=1000;b=0.005; y1=2*x; y2=cos(b*x); [haxes,hline1,hline2]=plotyy(x,y1,x,y2,'semilogy','plot'); axes(haxes(1)) ylabel('semilog plot'); axes(haxes(2)) ylabel('linear plot'); 实例5：单个轴窗口显示多个图形 function shili05 h0=figure('toolbar','none',... 'position',[200 150 450 250],... 'name','实例05'); t=0:pi/10:2*pi;

数值分析matlab外调函数

function [t,y]=chengmi(a,xinit,ep) v0=xinit; [tv,ti]=max(abs(v0)); lam0=v0(ti); u0=v0/lam0; flag=0; while (flag==0) v1=a*u0; [tv,ti]=max(abs(v1)); lam1=v1(ti); u0=v1/lam1; err=abs(lam1-lam0); if (err<=ep) flag=1; end lam0=lam1; end t=lam1; y=u0

-function [t,y]=fanchengmi(a,xinit,ep) A=a^(-1); v0=xinit; [tv,ti]=max(abs(v0)); lam0=v0(ti); u0=v0/lam0; flag=0; while (flag==0) v1=A*u0; [tv,ti]=max(abs(v1)); lam1=v1(ti); u0=v1/lam1; err=abs(lam1-lam0); if (err<=ep) flag=1; end lam0=lam1; end t=lam1 y=u0

3高斯消元法 function y=gaofun(A,b,n) m=ones(n,1); y=ones(n,1); for i=1:n for k=i+1:n m(k)=-A(k,i)/A(i,i); end for k=i+1:n for j=i:n A(k,j)=A(k,j)+m(k)*A(i,j); end b(k)=b(k)+m(k)*b(i); end end for k=n:(-1):1 s=0; for j=k+1:n s=s+A(k,j)*y(j); end

数值分析分段线性插值样条插值Runge函数Newton-Lagrange-Chebyshev插值多项式Runge现象matlab源程序代码

题目1：对R unge 函数R(x ) = 1 在区间[-1,1]作下列插值逼近，并和 x 2 R(x)的图像进行比较，并对结果进行分析。 = -1 + ih,h= 0.1,0 ≤ i≤ 20 绘出它的20 次Newton 插值(1)用等距节点x i 多项式的图像。 分别画出在[-1，1]区间，[-0.7，0.7]区间和[-0.5，0.5]区间上的 Newton 插值多项式和R unge 函数的图像

从图中可以看出： 1）在[-0.5，0.5]区间 N ewton 插值多项式和 R unge 函数的图像偏差较小，这说 明 Newton 插值多项式在此区间上可以较好的逼近 R unge 函数； 2）在边界（自变量 x =-1 和 x =1）附近，Newton 插值多项式和 R unge 函数的图像 偏差很大，Newton 插值多项式出现了剧烈的震荡。（Runge 现象） (2)用节点 x = cos(2i + 1 π）（, i = 0，1,2，? ? ? ，20），绘出它的 20 次 Lagrange i 42 插值多项式的图像。 画出在[-1，1]区间上的 L agrange 插值多项式和 R unge 函数的图像 从图中可以看出： 使用 C hebyshev 多项式零点构造的 L agrange 插值多项式和 R unge 函数的图 像偏差较小，没有出现 R unge 现象。

(3)用等距节点 x i 的图像。 = -1 + ih ,h = 0.1,0 ≤ i ≤ 20 绘出它的分段线性插值函数 画出在[-1，1]区间上分段线性插值函数和 R unge 函数的图像 从图中可以看出： 使用分段线性插值函数和 Runge 函数的图像偏差较小，也没有出现 Runge 现象，只在自变量 x =0 处有稍许偏差。 (4)用等距节点 x i 函数的图像。 = -1 + ih ,h = 0.1,0 ≤ i ≤ 20 绘出它的三次自然样条插值 画出在[-1，1]区间上三次自然样条插值函数和 R unge 函数的图像 从图中可以看出： 使用三次自然样条插值函数和 Runge 函数的图像偏差较小，也没有出现 Runge 现象。

数值分析作业MATLAB

1.用二分法解方程 x-lnx=2 在区间【2 ，4】内的根方法: 二分法 算法： f=inline('x-2-log(x)'); a=2;b=4;er=b-a; ya=f(a); er0=.00001; while er>er0 x0=.5*(a+b); y0=f(x0); if ya*y0<0 b=x0; else a=x0; ya=y0; end disp([a,b]);er=b-a;k=k+1; end 求解结果： >> answer1 3 4 3.0000 3.5000 3.0000 3.2500 3.1250 3.2500 3.1250 3.1875 3.1250 3.1563 3.1406 3.1563 3.1406 3.1484

3.1445 3.1484 3.1445 3.1465 3.1455 3.1465 3.1460 3.1465 3.1460 3.1462 3.1461 3.1462 3.1462 3.1462 3.1462 3.1462 3.1462 3.1462 3.1462 3.1462 最终结果为: 3.1462 2.试编写MATLAB 函数实现Newton 插值，要求能输出插值多项式。对函数 1 41)(2+=x x f 在区间[-5，5]上实现10次多项式插值。 Matlab 程序代码如下： %此函数实现y=1/(1+4*x^2)的n 次Newton 插值，n 由调用函数时指定 %函数输出为插值结果的系数向量（行向量）和插值多项式 算法： function [t y]=func5(n) x0=linspace(-5,5,n+1)'; y0=1./(1.+4.*x0.^2); b=zeros(1,n+1); for i=1:n+1 s=0; for j=1:i t=1; for k=1:i if k~=j

数值分析中求解线性方程组的MATLAB程序(6种)

数值分析中求解线性方程组的MATLAB程序（6种） 1.回溯法（系数矩阵为上三角） function X=uptrbk(A,B) %求解方程组，首先化为上三角，再调用函数求解 [N,N]=size(A); X=zeros(N,1); C=zeros(1,N+1); Aug=[A B]; for p=1:N-1 [Y,j]=max(abs(Aug(p:N,p))); C=Aug(p,:); Aug(p,:)=Aug(j+p-1,:); Aug(j+p-1,:)=C; if Aug(p,p)==0 'A was singular.No unique solution.' break; end for k=p+1:N m=Aug(k,p)/Aug(p,p); Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1); end end D=Aug; X=backsub(Aug(1:N,1:N),Aug(1:N,N+1)); 2.系数矩阵为下三角 function x=matrix_down(A,b) %求解系数矩阵是下三角的方程组 n=length(b); x=zeros(n,1); x(1)=b(1)/A(1,1); for k=2:1:n x(k)=(b(k)-A(k,1:k-1)*x(1:k-1))/A(k,k); end 3.普通系数矩阵（先化为上三角，在用回溯法） function X=uptrbk(A,B) %求解方程组，首先化为上三角，再调用函数求解 [N,N]=size(A); X=zeros(N,1); C=zeros(1,N+1); Aug=[A B]; for p=1:N-1

matlab 数值解

matlab 数值解 Matlab 数值解 Matlab 是一种强大的数学软件，它包含了很多数学工具箱，可以用于数值分析和求解数学问题。在本文中，我们将介绍Matlab 中的数值解方法，包括数值积分、数值微分、非线性方程求解和常微分方程的数值解法。 数值积分 数值积分是一种数学方法，用于求解函数的定积分。在Matlab 中，可以使用 quad 和 quadl 函数进行数值积分。其中，quad 函数用于计算一般积分，而 quadl 函数用于计算不定积分。 数值微分 数值微分是一种数学方法，用于计算函数的导数。在Matlab 中，可以使用diff 和gradient 函数进行数值微分。其中，diff 函数用于计算一维函数的导数，而 gradient 函数用于计算多维函数的梯度。 非线性方程求解 非线性方程是一种形式为 f(x)=0 的方程，其中 f(x) 是一个非线性函数。在 Matlab 中，可以使用 fzero 和 fsolve 函数进行非线性方程求解。其中，fzero 函数用于求解单变量非线性方程，而fsolve 函

数用于求解多变量非线性方程。 常微分方程的数值解法 常微分方程是一种形式为y'=f(t,y) 的方程，其中y 是未知函数，t 是自变量，f(t,y) 是已知函数。在Matlab 中，可以使用ode45 和ode23 函数进行常微分方程的数值解法。其中，ode45 函数是一种常用的数值解法，可以求解大部分常微分方程，而 ode23 函数则是一种高效的数值解法，适用于求解简单的常微分方程。 总结 在本文中，我们介绍了Matlab 中的数值解方法，包括数值积分、数值微分、非线性方程求解和常微分方程的数值解法。这些方法可以帮助我们快速、准确地求解数学问题，提高数学建模的效率和精度。

MATLAB语言及应用数值分析

MATLAB语言及应用数值分析 MATLAB是一种高级的计算机编程语言以及环境。它主要用于数值分 析和科学计算等领域，并且它还广泛应用于科学和工程领域。MATLAB提 供了丰富的数值计算函数和工具箱，能够轻松地进行数据处理、图形绘制、数值模拟以及算法开发等任务。 MATLAB的语法简洁明了，易于学习和理解。它具有强大的数值分析 能力，可以处理各种数学问题，如方程求解、数值积分、微分方程求解等。MATLAB的优势之一是它支持矩阵运算，这使得对线性代数和计算机图形 学等问题的解决变得更加高效。 MATLAB的应用领域广泛，可以应用于工程、科学、金融、生物医学、图像处理、机器学习等领域。在工程领域，MATLAB被广泛用于系统建模 与仿真、信号处理、控制系统设计等方面。在科学研究中，MATLAB可以 用于数据分析、实验数据处理、数据可视化等任务。在金融领域，MATLAB 可以进行风险管理、投资组合优化、期权定价等操作。在生物医学领域，MATLAB可以进行图像处理、信号处理、生物信号分析等工作。在机器学 习领域，MATLAB提供了各种强大的工具箱，可以进行模式识别、分类、 聚类等任务。 MATLAB还具有易于扩展、灵活性和跨平台等优点。用户可以通过编 写自定义函数和脚本来扩展MATLAB的功能。MATLAB还支持与其他编程语 言的交互，如C、C++、Java等。此外，MATLAB可以在Windows、Linux 和Mac等操作系统上运行，便于不同平台之间的代码的移植和共享。 总的来说，MATLAB是一种功能强大的编程语言和环境，广泛应用于 数值分析和科学计算等领域。它具有简洁的语法、丰富的数值计算函数和