北京科技大学数值计算方法实验代码

前言：

数值计算方法实验可以使用Matlab、C/C++、Python、Java等语言进行编程，考虑到同学期数学实验课程使用Matlab进行，建议提前熟悉Matlab编程（也效率更高）。

本文中各实验使用Matlab编程实现。

据观察，不同学年实验题目会有所变化，故本文将会介绍Matlab使用的一些基础知识，帮助理解。

（本文代码已添加注释）

Matlab使用入门

实验一：非线性方程求根

实验二：线性方程组求解

实验三：数值积分

Matlab使用入门

变量区分字母大小写 ; 变量名最多可包含63个字符（字母、数字和下划线），而且第一个字符必须是英文字母

ctrl + C 暂停

clc 清除命令行所有数据

clear 清除所有变量

clear + X 清除X变量

who 、whos 可以查看当前所有变量

while、if 循环、选择需要使用end 结束

disp() 输出打印

^ 乘方

log(x) lnx

1.0e-n 10的负n次方

diff() 求函数导数

如果语句后面不加分号（;），matlab会以交互式来执行程序，每执行一步，都会输出，输出的结果显示在命令行窗口。

使用分号，matlab会执行这个语句，并且继续执行，除非主动打印结果，否则不会显示到屏幕上。

函数定义：
在matlab里面输入（或页面左上角点击新建脚本并重命名）： edit fun.m （xxx.m）
在弹出的窗口输入以下内容并保存

function f=fun(x,y)

f=x^2+sin(x.*y)+2*y;

格式：function (自变量) function <因变量>=<函数名>(自变量)

脚本(.m)文件能保存、能运行

建议写一个定义main.m脚本文件用于主程序输出

矩阵输入：

A=[1,2,5 ; 8,4,5 ; 5,5,9]

以逗号(,) 为界限分列、分号(;)为界限分行

inv(X)对矩阵求逆。
如果A是非奇异方阵，则B/A = B*inv(A)，A\B = inv(A)*B。/表示右除，\表示左除。
注意：使用inv时，必须对象为方阵。

triu(a,k) 求上三角阵，a为原矩阵，k为与对角线的斜距离
tril(a,k) 求下三角阵，a为原矩阵，k为与对角线的斜距离
//k为正为上距离，k为负为下距离

plot：生成函数图像
x=linspace(0,1);

// linspace(a,b,n)生成从a到b之间线性分布的n个元素的数组，如n省略则默认为100。
y=exp(x)+3*x.^3-2*x.^2+log(x)-1;
plot(x,y)

实验一：非线性方程求根

共含六个脚本(.m)文件，实现二分法、简单迭代法、牛顿迭代法

f.m%原函数为f=exp(x)+3*x^3-2*x^2+log(x)-1;
function f=fun(x)
f=exp(x)+3*x^3-2*x^2+log(x)-1;
endg.m%牛顿迭代法 x1=x-f(x)/f`(x) 公式对应函数
%原函数的导数为f'(x)=exp(x)+9*x^2-4*x+1.0/x
function g = g(x)
g=x-fun(x)/(exp(x)+9*x^2-4*x+1.0/x);
end
%或使用如下公式(可以直接求导数)
%syms a;
%f = a - (fun(a)./diff(fun(a)));
%g = subs(f,x);%牛顿迭代公式main.m%采用两种以上不同的算法求解exp(x)+3*x^3-2*x^2+ln(x)-1=0在[0,1]范围内的一个根，要求两次迭代误差小于1.0e-4。
%计算时将main.m中代码贴入命令行窗口
x0=0;
x1=0.5;
x2=1;
eps=1.0e-4;%计算精度
first_fun(x0,x2,eps);%二分法迭代
second_fun(x1,eps);%简单迭代法迭代
third_fun(x2,eps);%牛顿迭代法迭代
%plot:生成函数图像
x=linspace(0,1);
y=exp(x)+3*x.^3-2*x.^2+log(x)-1;
plot(x,y)first_fun.m%二分法计算函数的零点
%函数输出根的近似值和迭代次数
function [temp,k]=first_fun(down,up,eps)
k=0;
temp=0;
while abs(up-down)>eps  %循环迭代，未达精度要求继续执行if k>20breakendtemp=(up+down)/2;y3=fun(temp);if y3==0up=temp;down=temp;elseif fun(down)*y3<0up=temp;elsedown=temp;endk=k+1;%test%disp('二分法迭代次数='),disp(k);%disp('近似解=');double(temp)
end
disp('二分法迭代次数='),disp(k);
disp('近似解=');double(temp)second_fun.m%简单迭代法计算函数的零点
%函数输出根的近似值和迭代次数
function [k,x1]=second_fun(x0,eps)
k=1;
x1=fun(x0);
%test
%disp('简单迭代法迭代次数='),disp(k);
%disp('近似解=');double(x1)
while(abs(x1-x0)>eps)if k>20breakend    x0=x1;x1=fun(x1);k = k +1;%test%disp('简单迭代法迭代次数='),disp(k);%disp('近似解=');double(x1)
end
disp('简单迭代法迭代次数='),disp(k);
disp('近似解=');third_fun.m%牛顿迭代法计算函数的零点
%函数输出根的近似值和迭代次数
function [k,x1]=third_fun(x0,eps)
k=1;
x1=g(x0);
%test
%disp('牛顿迭代法迭代次数='),disp(k);
%disp('近似解='),disp(x1);
while abs(x0-x1)>eps %eval(x0-x1)x0=x1;if k>20breakendx1 =g(x0);k = k+1;%test%disp('牛顿迭代法迭代次数='),disp(k);%disp('近似解='),disp(x1);
end
disp('牛顿迭代法迭代次数='),disp(k);
disp('近似解='),disp(x1);

实验二：线性方程组求解

共含五个脚本(.m)文件，实现高斯消元法、列主元消元法、雅可比迭代法、高斯-赛德尔迭代法

main.m%已知如下四元一次线性方程组,请选择至少两种算法编程求解它，并分析所得到结果的合理性。
a=[14,2,1,5;8,17,2,10;4,18,3,6;12,26,11,20];%原系数矩阵
b=[1;2;3;4];%常数列向量
x0=[0;0;0;0];
eps=1.0e-6;%精度
%disp('高斯消元法')
c=gauss(a,b);
%disp('列主元消元法')
d=lzy(a,b);
%disp('Jacobi迭代法')
[e,n1]=Jacobi(a,b,x0,eps);
%disp('高斯赛德尔迭代法')
[f,n2]=gauss_seidel(a,b,x0,eps);gauss.m%高斯消元法求解
function x=gauss(a,b)
%参量说明：A为系数矩阵；B为常数列向量；extend为增广矩阵
%将增广矩阵化为上三角，再回带求解x
%此方法较为常规，将extend(k,k)元素乘以-extend(i,k)/extend(k,k)加到第i行
%从1:n-1列，主对角元素的以下行，通过两层循环来遍历
extend=[a,b];
n=length(b);
ra=rank(a);
rz=rank(extend);
temp=rz-ra;
if temp>0disp('无解');return;
end
if ra==rzif ra==nx=zeros(n,1);for p=1:n-1for k=p+1:nm=extend(k,p)/extend(p,p);extend(k,p:n+1)=extend(k,p:n+1)-m*extend(p,p:n+1);endendb=extend(1:n,n+1);a=extend(1:n,1:n);x(n)=b(n)/a(n,n);for q=n-1:-1:1x(q)=(b(q)-sum(a(q,q+1:n)*x(q+1:n)))/a(q,q);endend
end
disp('高斯消元法解得:'),disp(x);lzy.m%列主元消元法求解
function x=lzy(a,b)
matrix=[a,b];
[row,col]=size(matrix);
x=zeros(row,1);%判定输入矩阵是否符合要求
if row~=col-1disp('请输入n*(n+1)维矩阵');
elsefor ii = 1:row-1%寻找主元max = ii;for jj = ii+1:rowif abs(matrix(jj,ii)) > abs(matrix(max,ii))max = jj;endendmatrix([ii,max],:) = matrix([max,ii],:);if matrix(ii,ii) == 0disp(['第',num2str(ii),'个主元素为零']);return;end%开始消元for jj = ii+1:rowmatrix(jj,:) = matrix(jj,:)-...matrix(jj,ii)/matrix(ii,ii)*matrix(ii,:);endend%消元完毕，开始回代if matrix(row,row)==0disp(['第',num2str(row),'个主元素为零']);return;endx(row)=matrix(row,col)/matrix(row,col-1);for ii=row-1:-1:1x(ii)=(matrix(ii,col)...-matrix(ii,1:row)*x)/matrix(ii,ii);end
end
disp('列主元消元法解得:'),disp(x);
endJacobi.m%Jacobi迭代法
function [y,n]=Jacobi(a,b,x0,eps)D=diag(diag(a));     %求对角矩阵
L=-tril(a,-1);       %求下三角矩阵
U=-triu(a,1);        %求上三角矩阵
B=D\(L+U);
f=D\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=eps %用2范数去逼近x0=y;y=B*x0+f;n=n+1;
end
disp('Jacobi迭代法解得:'),disp(y);
disp('迭代次数：'),disp(n);gauss_seidel.m%高斯赛德尔迭代法
function [y,n]=gauss_seidel(a,b,x0,eps)
D=diag(diag(a));   %求对角矩阵
L=-tril(a,-1);     %求下三角阵
U=-triu(a,1);      %求上三角阵
B=(D-L)\U;
f=(D-L)\b;
y=B*x0+f;
n=1;
while norm(y-x0)>=eps  %用2范数去逼近x0=y;y=B*x0+f;n=n+1;
end
disp('GaussSeidel迭代法解得:'),disp(y);
disp('迭代次数：'),disp(n);

实验三：数值积分

共含四个脚本(.m)文件，实现复合梯形公式、复合Simpson公式

fun.m%原函数
function f = fun(x)
f=(2*x-x^2+x^3+exp(x/2))^0.5;
endmain.m%选择一种或多种数值积分方法求解下列函数在区间[0,2]内的积分值，要求精度小于1.0e-5：
%f(x) =(2*x-x^2+x^3+exp(x/2))^0.5
a=0;
b=2;
e=1.0e-5;%精度
n=1;
Simpson(a,b,n,e);
tixing(a,b,n,e) ;
%plot:生成函数图像
%x=linspace(0,2);
%y=(2*x-x.^2+x.^3+exp(x/2)).^0.5;
%plot(x,y)tixing.mfunction [h,n] = tixing(a,b,n,e) %复合梯形公式
%区间为[a,b],n为初始区间切割次数,e为精度
%先求初始步长对应的公式值
h = (b-a)/n; %步长
T1= fun(a)+fun(b);
T1=T1*h/2;
%求变步长对应的公式值
n=n+1;
h = (b-a)/n;
x=a;
T2=fun(x)+fun(b);
for k=1:n-1       x = x+h;T2 = T2+ 2*fun(x);
end
T2=T2*h/2;
while abs(T2-T1)>e %循环迭代，未达精度要求继续执行T1=T2;x=a;T2=fun(x)+fun(b);n=n+1;h = (b-a)/n;for k=1:n-1x = x+h;T2 = T2+ 2*fun(x);endT2=T2*h/2;
end
disp('复合梯形公式迭代积分值:'),disp(T2);
disp('最终步长:'),disp(h);
disp('迭代次数:'),disp(n);Simpson.mfunction [h,n] = Simpson(a,b,n,e) %复合Simpon公式
%区间为[a,b],n为初始区间切割次数,e为精度
%先求初始步长对应的公式值
h = (b-a)/n; %步长
x = a;
S1= fun(x)-fun(b);
for k=1:n        x = x+h/2;S1 = S1+ 4*fun(x);x=x+h/2;S1=S1+ 2*fun(x);
end
S1=S1*h/6;
%求变步长对应的公式值
n=n+1;
h = (b-a)/n;
x=a;
S2=fun(x)-fun(b);
for k=1:n       x = x+h/2;S2 = S2+ 4*fun(x);x=x+h/2;S2 = S2+ 2*fun(x);
end
S2=S2*h/6;
while abs(S2-S1)>e  %循环迭代，未达精度要求继续执行S1=S2;x=a;S2=fun(x)-fun(b);n=n+1;h = (b-a)/n;for k=1:nx = x+h/2;S2 = S2+ 4*fun(x);x=x+h/2;S2 = S2+ 2*fun(x);endS2=S2*h/6;
end
disp('复合Simpon公式迭代积分值:'),disp(S2);
disp('最终步长:'),disp(h);
disp('迭代次数:'),disp(n);

北京科技大学数值计算方法实验代码相关推荐

北科大matlab第一次作业,北京科技大学应用计算方法作业与答案
<北京科技大学应用计算方法作业与答案>由会员分享,可在线阅读,更多相关<北京科技大学应用计算方法作业与答案(16页珍藏版)>请在人人文库网上搜索. 1.一.第一次作业(一)2- ...
北京科技大学计算机控制大作业,北京科技大学计算机控制系统实验报告
北京科技大学计算机控制系统实验报告 (23页) 本资源提供全文预览,点击全文预览即可全文预览,如果喜欢文档就下载吧,查找使用更方便哦! 17.90 积分计算机控制技术课程计算机控制技术课程实验报告 ...
北京科技大学工科物理实验大二上
前言本文由20级学生整理,包括实验目的和仪器.实验原理.实验步骤三个部分.主要是想节约一下大家手机拍照扫描.语音输入或手打的时间.(可能有些任课老师要求手写,那就爱莫能助了) 使用方法点击&quo ...
青岛科技大学c语言实验报告,青岛科技大学大学物理实验报告
青岛科技大学大学物理实验报告Tag内容描述: 1.北京科技大学实验报告磁场分布实验目的原理及实验步骤见预习报告实验数据附后及其处理 1 不同磁极头间隙内的磁场分布特点情形如图所示根据 ...
数值计算方法（停更笔记）
数值计算方法是通过计算机计算数学公式的一门科目,在以后的工作学习中,可能会常用这些方法.包括线性方程.非线性方程.数值积分.数值微分.微分方程的解.最小二乘法--在学习过程中就整理好这些方法的代码实现 ...
数值计算方法 | C语言实现几个数值计算方法（实验报告版）
目录写在前面实验一牛顿插值方法的实现实验二龙贝格求积算法的实现实验三高斯列主元消去法的实现实验四最小二乘方法的实现写在前面使用教材:<数值计算方法>黄云清等编著科学 ...
数值计算方法在计算机的应用,数值计算方法在计算机科学中的应用和误差序列实验推荐.doc...
数值计算方法在计算机科学中的应用和误差序列实验推荐数值计算方法在计算机科学中的应用和误差序列实验 [摘要]计算数学也叫做数值计算方法或数值分析.主要内容包括代数方程.线性代数方程组.微分方程的数值解 ...
数值计算方法上机c语言编程,数值计算方法上机实验报告.doc-资源下载在线文库www.lddoc.cn...
<数值计算方法>上机实验报告.doc 华北电力大学实验报告实验名称数值计算方法上机实验课程名称数值计算方法专业班级电力实 08 学生姓名李超然学号20080100 ...
matlab的数值求解实验报告,matlab计算方法实验报告5(数值积分)
计算方法实验报告(5) 学生姓名杨贤邦学号指导教师吴明芬实验时间2014.4.16地点综合实验大楼203 实验题目数值积分方法实验目的●利用复化梯形.辛普森公式和龙贝格数值积分公式计算定积分的近似 ...

北京科技大学数值计算方法实验代码

Matlab使用入门

实验一：非线性方程求根

实验二：线性方程组求解

实验三：数值积分

北京科技大学数值计算方法实验代码相关推荐

最新文章

热门文章

北京科技大学 数值计算方法实验代码

Matlab使用入门

实验一：非线性方程求根

实验二：线性方程组求解

实验三：数值积分

北京科技大学 数值计算方法实验代码相关推荐

最新文章

热门文章

北京科技大学数值计算方法实验代码

北京科技大学数值计算方法实验代码相关推荐