数值分析实验五
第六章 常微分方程数值解
1. 实验目的
1) 熟悉Matlab矩阵引用及高维矩阵操作 2) 用Matlab实现数值微分Richardson外推算法 3) 实现4阶Runge_Kutta常微分方程(组)数值解 4) 应用 Runge_Kutta方法解常微分方程(组)
2. 实验内容
(1) 数值微分Richardson外推算法
Richardson外推加速法递推公式: Tk4mk?1m?4m?1Tm?1?14m?1Tkm?1 算法实现:
function [val]=Jacbi_div(f, ,n,h) % Richardson外推加速法实现数值微分 %[val]=Jacbi_div(f, ,n,h)
% val 返回求导 (向量)函数的Jacbi矩阵 % f 求导 (向量)函数名
% 自变量(向量)值(不支持向量运算) % n 最大跌代次数(默认3) % h 初始步长 (默认0.1) if nargin<=3 h=0.1; end if nargin==2 n=3; end
r=size( ); m=r(2); b=size(feval(f, )); b=b(2); df=zeros(b,m); for t=1:m
Y= ; Y(t)=Y(t)+h; Z= ; Z (t)=Z (t)-h;
21
(k?1,2,3)
df(:,t)=((feval(f,Y)-feval(f,Z))./(2*h))'; end G(1,1,:,:)=df; for t=2:n h=h/2;
G=[G zeros(t-1,1,b,m);zeros(1,t,b,m)]; for k=1:m
Y= ; Y(k)=Y(k)+h; Z= ; Z (k)=Z (k)-h;
df(:,k)=((feval(f,Y)-feval(f,Z))./(2*h))'; end G(t,1,:,:)=df; for k=2:t r=k-1;
G(t,k,:,:)=((4^r).*G(t,r,:,:)-G(t-1,r,:,:))./(4^r-1); end % Richardson外推加速 end
val=reshape(G(n,n,:),[b,m]); %将数据整理成矩阵
例1: 求向量函数
y1?1?sin2x1?ex2 y2?x13?cos?x2?x1?
y3?x1?x2?x3在(0.9,2,1),(2.3,1.2,0.82)点的Jacbi矩阵 步骤:
1) 新建函数fun1.m function y=fun1(x)
y(1)=sqrt(1-sin(x(1)).^2)+exp(x(2)); y(2)=x(1).^3+cos(x(2)+x(1)); y(3)=x(1)+x(2)+x(3); 2)
>> [val]=Jacbi_div('fun1',[0.9 2 1],8,0.5) (2) 常微分方程(组)数值解--- Runge_Kutta方法
四阶 Runge_Kutta经典方法:
22
yn?1?yn?h?k1?2k2?2k?3k?;46k1?f?xn,yn?;
hh??k2?f?xn?,yn?k1?;22??hh??k3?f?xn?,yn?k2?;22??k4?f?xn?h,yn?hk1?;
算法实现:
function [T]=Runge_Kutta(f,x0,y0,h,n) % [T]=Runge_Kutta(f,x0,y0,h,n) % f 待解方程(组) % x0 初试自变量值 % y0 初试函数值 % h 步长 % n 步数 if nargin<5 n=100; end
r=size(y0); r=r(1); s=size(x0); s=s(1); r=r+s; T=zeros(r,n+1); T(:,1)=[y0;x0]; for t=2:n+1
k1=feval(f,T(1:r-1,t-1));
k2=feval(f,[k1*(h/2)+T(1:r-1,t-1);x0+h/2]); k3=feval(f,[k2*(h/2)+T(1:r-1,t-1);x0+h/2]); k4=feval(f,[k3*h+T(1:r-1,t-1);x0+h]); x0=x0+h;
T(:,t)=[T(1:r-1,t-1)+(k1+k2*2+k3*2+k4)*(h/6);x0]; end
23
数值分析实验六
第七章 非线性方程
1. 实验目的
1) 熟悉用Matlab实现迭代
2) 用Matlab实现二分,牛顿等求根算法 3) 应用算法求解方程根
2. 实验内容
(1) 二分法
二分法不断搜索有根区间,最终收缩为一点. 算法简单, 且收敛性总能得到保障. 实现如下:
function [val,n]=ErFen_root(f,x,dalt)
% [val,n]=ErFen_root(f,x,dalt) % f 要求根函数名 % x 初试有根区间 % dalt 精度,默认为 10 % val 返回所得根 % n 跌代次数 if nargin<3 dalt=1e-5; end
ab=x;i=0;x0=sum(ab)/2; f0=feval(f,ab(1)); f1=feval(f,x0); while abs(f1)>=dalt if f0*f1<0 ab(2)=x0; else ab(1)=x0; end
x0=sum(ab)/2; f0=feval(f,ab(1)); f1=feval(f,x0); i=i+1;
24
disp(sym([ab(1) x0 ab(2) f0 f1 feval(f,ab(2))])); % 显示迭代过程 if (ab(2)-x0)<2*eps
disp('May not be root!'); break; end end val=x0; if nargout==2 n=i; end
(2) Newton 切线法求非线性方程(组) Newton 切线法迭代公式:
? x1?x0?f0/f0
算法实现:
function [val,n]=Newton_root(f,a,daltf,daltx) % [val,n]=Newton_root(f,a,daltf,daltx) % val:输出解
% n:迭代次数(只有一个输出参数时不输出) % f:待解方程(组) % a:初始迭代值
% dlatf:函数精度控制(默认为1e-6) % daltx:解精度控制(默认为1e-5) tic %记时开始 if nargin<4 %设置默认值 daltx=1e-5; end if nargin<3 daltf=1e-6; end
flag=0; i=size(a);
if i(2)==1; %如果不为方程组,画图 flag=1; end
25