数值分析上机2new.docx
《数值分析上机2new.docx》由会员分享，可在线阅读，更多相关《数值分析上机2new.docx（18页珍藏版）》请在冰豆网上搜索。
数值分析上机2new
实验报告
（二）
课程名称《数值分析》
实验项目常微分方程的差分法、方程求根迭代法和方程组求解的迭代法
实验环境PC机、Matlab、C，C++
分组张申、王星
班级/学号/姓名电信10032010010552张申
指导教师王灯山
实验日期2012-12-04
成绩，
一、实验名称：
常微分方程的差分法、方程求根迭代法和方程组求解的迭代法
二、实验目的：
1、掌握求解一阶常微分方程初值问题的欧拉法、改进的欧拉法、龙格-库塔法；掌握方程求根的二分法、牛顿下山法、弦截法；分段插值方法、会用复化求积法；掌握方程组求解的Jacobi法等
2、理解各种插值格式的精度的理论分析；
3、能正确地将理论分析算法过程用Matlab或C相关的程序设计语言来实现；
4、对程序运行过程中出现的错误能自行调试修正；
5、能从数值结果得到关于精度的验证.
三、实验内容：
1.本题两人一组。
利用改进的欧拉法、变形欧拉法、三阶龙格-库塔格式、四阶龙格-库塔格式求解如下常微分方程（要求：
每人至少用两种方法求解，且其中一个必须是四阶龙格-库塔格式），并画出图形。
（2）
张申+王星%改y
（1）=1x的取值范围[1，2]
（一）三阶龙格库塔法
（1）源程序
clearall
a=1;%定义区间a，b
b=2;
N=10;%定义所分份数10份
ya=1;
h=（b-a）/N;%计算每一份所得h
y
（1）=ya;
x=a:
h:
b;%令x等于a与b分成h份
fori=1:
N%利用三阶龙格库塔法公式
K1=f2（x（i）,y（i））;
K2=f2（x（i）+0.5*h,y（i）+h/2*K1）;
K3=f2（x（i）+h,y（i）+h*（（-K1）+2*K2））;
y（i+1）=y（i）+h/6*（K1+4*K2+K3）;
end
yy=x;%得出精确解
plot（x,y,'rs',x,yy）%画出xy以及yy
[y,yy]
-----------------------------------------------------------
functionz=f2（x,y）%所计算的式子
z=x/y;
（2）结果图：
（二）四阶龙格库塔
（1）源程序（c语言）
#include"stdio.h"
#include"conio.h"
floatfunc（floatx,floaty）//定义浮点型xy
{
return（x/y）;//返回所算例题x/y
}
floatrunge_kutta（floatx0,floatxn,floaty0,intn）
{
floatx,y,y1,y2,h,xh;//定义各个变量
floatd1,d2,d3,d4;//定义各个变量
inti;//定义整形量
x=x0;//定义区间
y=y0;
h=（xn-x0）/n;//计算区间被分成n分
for（i=1;i<=n;i++）//利用for语句连续计算
{
xh=x+h/2;//公式
d1=func（x,y）;
d2=func（xh,y+h*d1/2.0）;
d3=func（xh,y+h*d2/2.0）;
d4=func（xh,y+h*d3）;
y=y+h*（d1+2*d2+2*d3+d4）/6.0;
x=x0+i*h;
}
return（y）;//返回y
}
voidmain（）//主程序
{
floatx0,xn,y0,e;
intn;
printf（"\ninputn:
\n"）;//输入n（所分份数）
scanf（"%d",&n）;//显示n
printf（"inputx0,xn:
\n"）;//输入区间
scanf（"%f%f",&x0,&xn）;//显示区间
printf（"inputy0:
\n"）;//输入y0
scanf（"%f",&y0）;//显示y0
e=runge_kutta（x0,xn,y0,n）;//返回子函数runge_kutta并计算
printf（"y（%f）=%6.6f",y0,e）;//输出y，
}
2.本题两人一组。
利用二分法、牛顿法、牛顿下山法和弦截法求解如下方程（要求：
每人至少用两种方法求解，且两人所用方法不能完全重复）。
要求精度为10^（-5）.
（2）
在x=0.5附近的一个根。
张申+王星
方法一二分法
源程序：
------------------------------------------------------------------------------
functionx=agui_bisect（fname,a,b,e）
%fname为函数名，a,b为区间端点，e为精度
fa=feval（fname,a）;%把a端点代入函数，求的fa
fb=feval（fname,b）;%把b端点代入函数，求的fb
iffa*fb>0error（'两端函数值为同号'）;end%如果fa*fb>0，则输出两端函数值为同号
k=0
x=（a+b）/2
while（b-a）>（2*e）%循环条件的限制
fx=feval（fname,x）;%把x代入代入函数，求的fx
iffa*fx<0%如果fa与fx同号，则把x赋给b，把fx赋给fb
b=x;
fb=fx;
else%如果fa与fx异号，则把x赋给a，把fx赋给fa
a=x;
fa=fx;
end
k=k+1%计算二分了多少次
x=（a+b）/2%当满足了一定精度后，跳出循环，每次二分，都得新的区间断点a和b，则近似解为x=（a+b）/2
end
1、在MATLAB命令窗口求解方程f（x）=
，即输入如下
>>fun=inline（'exp（-x）-x'）%所求算式
在x=0.5
>>x=agui_bisect（fun,0,1,10^-5）%0,1，精度10^-5
2、得到计算结果，且计算结果为（截图顺序从左向右）
结论：
二分了16次后达到精度x=0.5671
方法二牛顿法
matlab源程序：
第一步：
在MATLAB7.0软件，建立一个实现牛顿迭代法的MATLAB函数文件=agui_newton.m如下：
functionx=agui_newton（fname,dfname,x0,e）
%fname为函数名dfname的函数fname的导数，x0为迭代初值
%e为精度，N为最大迭代次数（默认为100）
N=100;
x=x0;%把x0赋给x，再算x+2*e赋给x0
x0=x+2*e;
k=0;
whileabs（x0-x）>e&kx0-x的绝对值大于某一精度，和迭代次数小于N
k=k+1%显示迭代的第几次
x0=x;
x=x0-feval（fname,x0）/feval（dfname,x0）;%牛顿迭代公式
disp（x）%显示x
end
ifk==Nwarning（'已达到最大迭代次数'）;end%如果K=N则输出已达到最大迭代次数
第二步：
在MATLAB命令窗口求解方程f（x）=e^x+10x-2=0，即输入如下
>>fun=inline（'exp（-x）-x'）%输入所计算式子
>>dfun=inline（'-exp（-x）-1'）%输入所计算式子导数
>>x=agui_newton（fun,dfun,0,10^-5）
第三步：
得出结果，且结果为
结论：
x=0.5671
3.本题5人一组。
利用Jacobi迭代法、Gauss-Seidel迭代或者超松弛迭代法求解如下线性方程。
要求精度为10^（-6）.
（3）.
张申+王星+罗伟+孙昱
Jacobi迭代法源程序：
function[x,k,index]=Jacobi（A,b,ep,it_max）
%求解线性方程组的Jacobi迭代法，其中
%A---方程组的系数矩阵
%b---方程组的右端项
%ep---精度要求。
省缺为1e-5
%it_max---最大迭代次数，省缺为100
%x---方程组的解
%k---迭代次数
%index---index=1表示迭代收敛到指定要求；
%index=0表示迭代失败
ifnargin<4it_max=100;end
ifnargin<3ep=1e-5;end
n=length（A）;k=0;
x=zeros（n,1）;y=zeros（n,1）;index=1;
while1
fori=1:
n
y（i）=b（i）;
forj=1:
n%for循环
ifj~=i
y（i）=y（i）-A（i,j）*x（j）;
end
end
ifabs（A（i,i））<1e-10|k==it_max
index=0;return;
end
y（i）=y（i）/A（i,i）;
end
ifnorm（y-x,inf）break;
end
x=y;k=k+1;
end
实验分析：
x1=x2=x3=-0.5
迭代次数为18次
index=1表示迭代收敛到指定要求
四、实验心得：
首先要向王老师道歉，曾经对这么课程的松懈，以至于错过了几次课程。
数值分析这门课程我并不太感兴趣，繁琐的数学问题非常枯燥。
但是数值分析这门课程的上机部分十分有趣，我学到了matlab的很多基础知识，当然也感谢小梁同学的帮助。
同时我也尝试着用c语言编写了几个程序。
今后再面对复杂繁琐的数学问题我便可以用电脑实现。
在今后的道路上我会更加重视数值分析以及matlab。
升级会员 