当前位置：文档之家 › Matlab(R2009a版)-第8讲 计算方法的MATLAB实现

Matlab(R2009a版)-第8讲 计算方法的MATLAB实现

  • 格式：pps
  • 大小：2.18 MB
  • 文档页数：169

MATLAB2009从入门到精通2013-7-222课程主要内容•第1章MATLAB简介•第2章数值运算•第3章单元数组和结构•第4章字符串•第5章符号运算•第6章MATLAB绘图基础•第7章程序设计•第8章计算方法的MATLAB实现•第9章优化设计•第10章SIMULINK仿真初探2013-7-223第8章计算方法的MATLAB实现•随着计算机的迅速发展与广泛运用，在众多的领域，科学计算方法的应用越来越广泛了，而MATLAB在进行科学计算方面有着无与伦比的优势。本章介绍MATLAB在计算方法中的运用。2013-7-2248.1 方程求根•roots见多项式求根；roots(多项式系数矩阵)•fzero可求解非线性方程，它的格式为：•fzero(‘function’,x0)•其中function为求解的方程，x0为估计的根，x0可为标量或长度为2的向量，为向量时函数的两端的值应该符号相反，此时求区间上的解。只能求解x0附近的一个解。即使在某个区间内有多个解，但是区间端点符号相同的话仍然出错。2013-7-225程序实例•>> fzero('x^3-3*x-1',2)•ans =•1.8794•>> fzero('x^3-3*x-1',[1,4])•ans =•1.8794•>> fzero('x^3-3*x-1',[2,4])•??? Error using ==> fzero•The function values at the interval endpoints must differ in sign.2013-7-226程序实例•>> fzero('x^2-3*x+2',0)•ans =•1.0000•>> fzero('x^2-3*x+2',3)•ans =•2.0000•>> fzero('x^2-3*x+2',[0,4])•??? Error using ==> fzero•The function values at the interval endpoints must differ in sign.2013-7-2278.2 线性方程组数值解法•线性方程组的求解不仅在工程技术领域涉及到，而且在其他的许多领域也经常碰到，因此这是一个应用相当广泛的课题。•关于线性方程组的数值解法一般分为两类：一类是直接法，就是在没有舍入误差的情况下，通过有限步四则运算求得方程组准确解的方法。另一类是迭代法，就是先给定一个解的初始值，然后按一定的法则逐步求出解的各个更准确的近似值的方法。2013-7-2288.2.1 直接解法•关于线性方程组的直接解法有许多种，比如Gauss消去法、列主元消去法、平方根法等。而在MATLAB中，线性方程组的直接解法只需用符号“/”或“\”就解决问题。还可以使用逆阵函数来求解：x=inv(A)*B。2013-7-229程序实例•求解下列方程组615318153312321321321xxxxxxxxx2013-7-2210程序实例•>> a=[12 -3 3;-18 3 -1;1 1 1];•>> b=[15;-15;6];•>> x1=a\b•x1 =•1.0000•2.0000•3.0000•>> x2=inv(a)*b•x2 =•1•2•32013-7-22118.2.2 线性方程组求解中的变换•上三角变换•U=triu(x)返回矩阵x的上三角部分；•U=triu(x,k)返回第k条对角线以上部分的元素。2013-7-2212程序实例•>> a=[12 -3 3;-18 3 -1;1 1 1];•>> triu(a)•ans =•12 -3 3•0 3 -1•0 0 12013-7-2213•>> triu(a,1)•ans =•0 -3 3•0 0 -1•0 0 0•>> triu(a,-1)•ans =•12 -3 3•-18 3 -1•0 1 1程序实例2013-7-2214•下三角变换•U=tril(x)返回矩阵x的下三角部分；•U=tril(x,k)返回第k条对角线以上下部分的元素。2013-7-2215程序实例•>> a=[12 -3 3;-18 3 -1;1 1 1];•>> tril(a)•ans =•12 0 0•-18 3 0•1 1 12013-7-2216程序实例•>> tril(a,1)•ans =•12 -3 0•-18 3 -1•1 1 1•>> tril(a,-1)•ans =•0 0 0•-18 0 0•1 1 02013-7-2217•对角变换•U=diag(x)返回矩阵x主对角线上的元素，返回结果是一列向量形式；•U=diag(x,k)返回第k条对角线上的元素值。•当x为向量时生成矩阵。2013-7-2218程序实例•>> a=[12 -3 3;-18 3 -1;1 1 1];•>> diag(a)•ans =•12•3•12013-7-2219程序实例•>> a=[12 -3 3;-18 3 -1;1 1 1];•>> diag(a,1)•ans =•-3•-1•>> diag(a,-1)•ans =•-18•12013-7-22208.2.3 迭代解法•迭代解法非常适用于求解大型稀疏系数矩阵的方程组，在线性方程组常用的迭代解法主要有Jacobi迭代法、Gauss-Seidel迭代法。2013-7-2221Jacobi迭代法2013-7-2222Jacobi.m•function s=jacobi(a,b,x0,eps)•%jacobi迭代法皆线性方程组•%a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值•if nargin==3•eps=1.0e-6;•elseif nargin<3•error•return•end2013-7-2223•D=diag(diag(a));%求出对角矩阵•D=inv(D);%求出对角矩阵的逆•L=tril(a,-1);%求出严格下三角矩阵•U=triu(a,1);%求出严格上三角矩阵•B=-D*(L+U);•f=D*b;•s=B*x0+f;•while norm(s-x0)>=eps•x0=s;•s=B*x0+f;•end•return2013-7-2224程序实例•用上面编写的jacobi函数求解下列方程组1052151023210321321321xxxxxxxxx2013-7-2225程序实例•>> a=[10 -2 -1;-2 10 -1;-1 -2 5];•>> b=[3 15 10]';•>> x=jacobi(a,b,[0 0 0]',eps)•x =•1•2•32013-7-2226Gauss-Saidel迭代法2013-7-2227gauss.m•function s=gauss(a,b,x0,eps)•%gauss-seidel迭代法皆线性方程组•%a为系数矩阵，b为方程组ax=b中的右边的矩阵b，x0为迭代初值•if nargin==3•eps=1.0e-6;•elseif nargin<3•error•return•end2013-7-2228•L=tril(a,-1);%求出严格下三角矩阵•D=diag(diag(a));%求出对角矩阵•U=triu(a,1);%求出严格上三角矩阵•C=inv(D+L);•B=-C*U;•f=C*b;•s=B*x0+f;•while norm(s-x0)>=eps•x0=s;•s=B*x0+f;•end•return2013-7-2229程序实例•用上面编写的gauss函数求解下列方程组105210226210321321321xxxxxxxxx2013-7-2230程序实例•>> a=[10 -2 -1;-2 2 -1;-1 -2 5];•>> b=[6;10;10];•>> x=gauss(a,b,[0 0 0]',eps)•x =•4•13•82013-7-22318.3 非线性方程组数值解法•与线性方程组的求解一样，非线性方程组的求解也是应用很广泛的课题。一般情况，非线性方程组的数据值解法主要采用迭代法来求解。比较常用的迭代法主要有不动点迭代法、Newton迭代法、拟Newton迭代法等几种方法。2013-7-2232不动点迭代法2013-7-2233staticiterate.m•function s=staticiterate(x,eps)•%不动点迭代法解非线性方程组，x为迭代初值，eps为允许误差•if nargin==1•eps=1.0e-6;•elseif nargin<1•error•return•end2013-7-2234•xx=fx(x);%第一次迭代•while norm(xx-x)>=eps•x=xx;•xx=fx(x);•end•s=xx;•return2013-7-2235程序实例•用不动点迭代法求解下面的方程组081008102122122121xxxxxxx2013-7-2236fx.m•首先编写上述非线性方程组的M文件fx.m•function y=fx(x)•y(1)=0.1*(x(1)*x(1)+x(2)*x(2)+8);•y(2)=0.1*(x(1)*x(2)*x(2)+x(1)+8);•y=[y(1) y(2)];2013-7-2237程序实例•>> x=staticiterate([0 0])•x =•1.0000 1.00002013-7-2238Newton迭代法2013-7-2239newtoniterate.m•function s=newtoniterate(x,eps)•%newton迭代法解非线性方程组，x为迭代初值，eps为允许误差•if nargin==1•eps=1.0e-6;•elseif nargin<1•error•return•end2013-7-2240•x1=fx1(x);%非线性方程组•x2=-dfx1(x);%非线性方程组导数•x3=inv(x2);•x0=x3*x1';•while norm(x0)>=eps•x=x0'+x;•x1=fx1(x);•x2=-dfx1(x);•x3=inv(x2);•x0=x3*x1';•end•s=x0'+x;•return2013-7-2241程序实例•用上面编写的newton迭代函数求解下列方程组081008102122122121xxxxxxx2013-7-2242fx1.m和dfx1.m•首先，编写上述非线性方程组的M文件fx1.m•function y=fx1(x)•y(1)=x(1)*x(1)-10*x(1)+x(2)*x(2)+8;•y(2)=x(1)*x(2)*x(2)-10*x(2)+x(1)+8;•y=[y(1) y(2)];•然后，编写上述非线性方程组导数的M文件dfx1.m•function y=dfx1(x)•y(1)=2*x(1)-10;•y(2)=2*x(2);•y(3)=x(2)*x(2)+1;•y(4)=2*x(1)*x(2)-10;•y=[y(1) y(2);y(3) y(4)];2013-7-2243程序实例•>> x=newtoniterate([0 0])•x =•1 12013-7-22448.4 插值与拟合•在生产实践及科学实验中，插值与拟合的应用非常广泛。下面，就对如何用MATLAB来处理插值与拟合作一介绍。2013-7-22458.4.1 一维插值•yi=interp1(x,y,xi)返回在插值向量xi处的函数向量yi，它是根据向量x和y插值而来。如y是矩阵，则对y每一列进行插值，如xi中元素不在x内，返回NaN。•yi=interp1(y,xi)省略x，表示x=1:N，此时N为向量y的长度或为矩阵y的行数。•yi=interp1(x,y,xi,’method’)表示用method指定的插值方法进行插值。

合集下载

matlab计算涡度的函数_流函数涡量法的二维方腔流数值模拟matlab编程.doc

matlab计算涡度的函数_流函数涡量法的二维方腔流数值模拟matlab编程.doc

matlab计算涡度的函数_流函数涡量法的⼆维⽅腔流数值模拟

matlab编程.doc

流函数- 涡量法的⼆维⽅腔流数值模拟

基本⽅程：

在直⾓坐标系下，不可压⾮定常流体所满⾜的流函数涡量形式的N-S⽅程为

其中

为雷诺数

差分格式：

采⽤FTCS格式有：

对于本问题，将⽅腔四边同时分为等分，则有

故

在直⾓坐标系下，不可压定常流体所满⾜的流函数涡量形式的N-S⽅程为

其中

为雷诺数

差分格式：

采⽤FTCS格式有：

对于本问题，将⽅腔四边同时分为等分，则有，则有即

边界条件：

在腔体的两侧和顶边，

(第⼆式由泰勒级数展开得到)

在底边

(第⼆式由泰勒级数展开得到)

其中代表边界，代表与边界相邻的节点。

⽽

即

Matlab程序为：

不可压⾮定常流体

clear;

%参数设置

Re=10; %雷诺数取10,100，500，1000

L=1; %空⽳⼏何尺⼨

n=100;

dh=L/n;%delta hdt=1e-4; %时间步长

psi=zeros(n+1,n+1);

xi=zeros(n+1,n+1);

rho=1;

for k=1:1000000

err=0;

%边界条件

for i=2:n

xi(i,1)=-2*(psi(i,2)-psi(i,1))/dh^2;

xi(i,n+1)=-2*(psi(i,n)-psi(i,n+1))/dh^2;

end

for j=2:n

xi(1,j)=-2*(psi(2,j)-psi(1,j)+dh)/dh^2;

xi(n+1,j)=-2*(psi(n,j)-psi(n+1,j))/dh^2;

end

%控制⽅程

for i=2:n

for j=2:n

u(i,j)=(psi(i,j+1)-psi(i,j-1))/(2*dh);

v(i,j)=-((psi(i+1,j)-psi(i-1,j))/(2*dh));

err1=(psi(i+1,j)+psi(i-1,j)+psi(i,j+1)+psi(i,j-1)+xi(i,j)*dh^2)/4-psi(i,j);

MATLAB电子教案

MATLAB电子教案

第 1 页/共 63页 第1章 MATLAB简介

1.1 概述

MATLAB是MATrix LABoratory（矩阵实验室）的缩写，由美国The Math Works公司于1984年推出的一种科学与工程计算语言。主要特点：

一、简单易学。

二、代码短小高效。

三、功能丰富，可扩展性强。

四、强大的图形表达功能。

五、强有力的系统仿真功能。

1.2 桌面启动

启动MATLAB桌面主要采用以下两种方法：

一、在Windows桌面上，双击MATLAB的快捷方式图标。采用这种方式打开的MATLAB桌面以matlab*\work为当前目录。

二、双击matlab*\bin\win32文件夹中的MATLAB.exe。采用这种方式打开的MATLAB桌面以matlab71为当前目录。

两者区别：当前目录不同。

注：*为MATLAB的软件版本号

1.3 通用操作界面简介

一、命令窗口

缺省情况下，位于桌面右侧，是用户与MATLAB进行人机对话的主要环境。在该窗口内，可输入各种由MATLAB运行的命令、函数、表达式，显示除图形外的所有运算结果。

二、命令历史窗口

缺省情况下，位于桌面左下方的前台，该窗口记录并显示每次开启MATLAB的时间及所有MATLAB运行过的命令、函数及表达式等，允许用户对它们进行选择复制、重运行及产生M文件。

三、当前目录浏览器

缺省情况下，位于MATLAB桌面左上方的前台。在该浏览器中，可以进行当前目录的设置，展示相应目录上的.m及.mdl等文件，复制、编辑和运行M文件以及装载MAT数据文件等。

四、工作空间浏览器

缺省情况下，位于MATLAB桌面左上方的后台，该窗口列出了MATLAB工作空间中所有数据的变量信息，包括变量名、大小、字节数等。在该窗口中，可以对变量进行观察、编辑、提取及保存。

五、数组编辑器

缺省情况下，不随操作界面的出现而启动，只有在工作空间浏览器中对变量进行操作时才启动。

matlab lsqcurvefit 拟合度r2计算函数

matlab lsqcurvefit 拟合度r2计算函数

在MATLAB中，`lsqcurvefit` 用于非线性最小二乘拟合，它不直接提供R^2（R方，也叫决定系数）的计算。R^2是用于衡量模型拟合优度的一个指标，值越接近1表示模型拟合越好。

你可以通过以下步骤计算R^2：

1. 使用`lsqcurvefit`得到拟合的参数。

2. 利用这些参数生成预测值。

3. 计算实际值和预测值之间的均方误差（MSE）。

4. 计算R^2。R^2等于1减去MSE除以总体的方差。

以下是一个示例代码：

```matlab

假设有以下数据

x = linspace(0, 10, 100); x数据

y = sin(x) + randn(size(x)); y数据，其中sin(x)是信号，randn(size(x)生成随机噪声

定义要拟合的函数形式，例如这里是一个简单的线性函数

fun = @(p,x) p(1)*x + p(2);

为拟合的参数初始化一个初始猜测

p0 = [0, 1];

使用lsqcurvefit进行拟合

p = lsqcurvefit(fun, p0, x, y);

生成预测值

y_pred = fun(p, x);

计算MSE

mse = mean((y - y_pred).^2);

计算R^2

r2 = 1 - mse / var(y);

输出R^2

fprintf('R^2: .2f\n', r2);

```

注意：在这个例子中，我假设了你正在拟合一个线性模型。如果你正在拟合一个非线性模型，那么你可能需要为`fun`定义一个更复杂的函数形式，并且可能需要更复杂的参数初始化策略。

数值计算方法matlab 第一章 误差分析

数值计算方法matlab 第一章 误差分析

数值计算方法

1 第一章作业

1.对一个数求和100000次。对数1以单精度方式求和，对数0.00001分别以单精度和双精度方式求和。

问题分析：单精度方式使用函数single()，双精度求和为matlab自动调整，不需要特别说明。

程序编写如下：

运行结果：

实验结果分析：

不难看出，对于1进行单精度求和得到的结果和期望值一致，但是对0.00001进行单精度求和的结果却存在误差，对0.00001进行双进度求和，误差得到减小。这是由于量化误差造成的，0.00001在计算机中并不能准确表示，只能对其进行量化处理，得到一个和真值有一点区别的量化值，小量计算中可以忽略，但在计算了100000后误差积累，导致了最后的结果误差较大。双精度的情况下，该误差小得多。 数值计算方法

2 2.如果||1x，23111xxxx…

当x=0.1时，从111x开始，然后每次加入一项来分别计算1/0.9。在每加入一个新项后，计算近似百分比相对误差，直到近似误差估计值的绝对值小于与五位有效数字一致的误差准则时停止计算。

问题分析：

本例中，要保证5位有效数字，因此容限误差为：

256s(0.510)%510

近似百分比误差为：

-100%a当前近似值前一近似值当前近似值

真误差为：

-100%真值近似值真值

跳出循环的标准为：a|s

程序编写如下：

数值计算方法

3 运行结果如下：

实验结果分析：

实验结果表明，当计算到第6次时，近似误差就已经小于了容限值，循环结束。随着添加多的项数，实际误差和近似误差都减小了，说明了计算精度在逐步提高。我们可以通过改变s的值来调节所需要的计算精度。