⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 m文件列表.txt

📁 (有源代码)数值分析作业,本文主要包括两个部分,第一部分是常微分方程(ODE)的三个实验题,第二部分是有关的拓展讨论,包括高阶常微分的求解和边值问题的求解(BVP).文中的算法和算例都是基于Matla
💻 TXT
字号:
 Volume in drive E has no label.
 Volume Serial Number is 2841-B3A8

 Directory of E:\matlabtemp

2008-05-01  09:42             1,078 adams_wjl.m
2008-05-03  15:58             1,267 adams_wjl1.m
2008-05-04  18:27               662 bvptest.m
2008-04-30  15:22                46 cd_temp.m
2008-04-30  15:23                50 cd_work.m
2008-05-04  12:48               181 dampmass.m
2008-05-04  12:47               184 dampmass1.m
2008-05-04  12:50               168 dampmass_test.m
2008-05-03  21:18               298 digui.m
2008-05-03  21:14               103 digui_r.m
2008-05-01  22:52               112 exact51.m
2008-05-01  11:27               204 exact52.m
2008-05-03  15:13               204 exact53.m
2008-04-30  17:01               158 funtest1.m
2008-05-03  15:58               102 funtest2.m
2008-05-03  15:58               478 funtest3.m
2008-05-01  22:54            23,968 history080501.m
2008-05-03  15:58               136 na52test.m
2008-05-03  22:47             2,752 na53test.m
2008-04-30  15:59             1,083 nm_ch_0501.m
2008-04-15  21:42               331 ODEfun5.m
2008-05-04  14:35                88 odesol_w.m
2008-04-30  16:55               720 RK1.m
2008-05-04  16:36             2,071 rk51test.m
2008-04-30  18:15               174 rk_ex_5_1.m
2008-04-30  18:10               641 rk_wjl.m
2008-05-04  19:29               140 vanderpolequ.m
2008-05-04  10:36             2,319 vdpdemo_wjl.m
2008-05-04  16:00               596 vdptest.m
2008-04-15  22:00               608 wjlvdp.m
              30 File(s)         40,922 bytes
               0 Dir(s)     366,952,448 bytes free


a文件列表:
(只列出了核心文件,其它一些零散的文件未列出)
adams_wjl1.m   Adams方法求解一阶常微分方程
bvptest.m      边界值问题求解例子
exact51.m      5.1问题的精确解函数
exact53.m      5.3 问题的精确解函数
funtest1.m     5.1问题所用到的函数
funtest2.m     5.2问题所用到的函数
funtest3.m     5.3问题所用到的函数
na52test.m     数值分析5.2的求解例子
na53test.m     数值分析5.3的求解例子
rk51test.m     数值分析5.1的求解例子
vanderpolequ.m VDP方程的函数
vdptest.m      VDP求解的例子

adams_wjl1.m   Adams方法求解一阶常微分方程
function [x,y]=adams_wjl1(odefile,xi,xf,yi,h,varargin)

% odefile is the filename of ode , xi<x <xf,  yi is the initial value, h is
% the step , varargin  is used by odefile

%
%  See also ODE23, ODE45, ODE113, ODE15S, ODE23S,
%           

x = [xi:h:xf];		% Vector of x values
if x(end) ~= xf
   x(end+1) = xf;
end

d = diff(x);
% if you don't have this sentence,change d(i) to h is OK

% Starting values
% tic;

% [a,b] = rk_wjl(odefile,x(1),x(4),h,yi,varargin{:});
[a,b] = rk_wjl(odefile,x(1),x(4),yi,h,varargin{:});
%do not write yi , h, in the wrong direction
y(1:4) = b;
for i = 1:4
   f(i) = feval(odefile,x(i),y(:,i),varargin{:});
end

ynp=0;
ync=0;



% Solution
for i = 4:length(x)-1
   % Predictor
   yn1p = y(i) + d(i)/24 * (55*f(i) - 59*f(i-1) ...
      + 37*f(i-2) - 9*f(i-3));
   yn1pm=yn1p +251/270*(ync-ynp);
   
   %f(i+1) = feval(odefile,x(i+1),y(i+1),varargin{:});
    f(i+1) = feval(odefile,x(i+1),yn1pm,varargin{:});

   % Corrector
   yn1c = y(i) + d(i)/24 * (9*f(i+1) + 19*f(i) ...
      - 5*f(i-1) + f(i-2));
   y(i+1) =yn1c -19/270*(yn1c - yn1p);
   %f(i+1) = feval(odefile,x(i+1),y(i+1),varargin{:});
   ync=yn1c;
   ynp=yn1p;
   
end
% time=toc;
% you use time out the programme

bvptest.m      边界值问题求解例子
% function bvptest
solinit = bvpinit([0 1 2 3 4],[100; 0]);
sol = bvp4c(@twoode, @twobc, solinit);

xint = linspace(0, 4, 50);
yint = deval(sol, xint);
plot(xint, yint(1,:),'b');
xlabel('x')
ylabel('solution y')
% hold on

%%
% This particular boundary value problem has exactly two solutions.  The other
% solution is obtained for an initial guess of 
% 
%     y(x) = -1, y'(x) = 0 
% 
% and plotted as before.

hold on
clear
solinit = bvpinit([0 1 2 3 4],[-100; 0]);
sol = bvp4c(@twoode,@twobc,solinit);

xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');

legend('solinit +100','solinit -100')
title('BVP')
exact51.m      5.1问题的精确解函数
function  y=exact51(x,a)
% the exact solution for  chapter 5 problem1 ,
% y = exp(a*x) +x;

y = exp(a*x) +x;

exact53.m      5.3 问题的精确解函数
function  y=exact53(x)
% the exact solution for  chapter 5 problem1 ,
% s= dsolve('Dy=-y +cos(2*x) - 2*sin(2*x) + 2*x*exp(-x)','y(0)=1','x')
% s =
% cos(2*x)+exp(-x)*x^2

y = cos(2*x)+exp(-x).*x.^2;

funtest1.m     5.1问题所用到的函数
function yy=funtest1(x,y,a)
% funtest1 
%yy=1 + (y - x) .* a;
%to call it  use


%feval('funtest1',1,2,0.5)
% ans =
%     1.5000
yy=1 + (y - x) .* a;

funtest2.m     5.2问题所用到的函数
function f=funtest2(x,y)
f = -y + cos(2*x) - 2*sin(2*x)...
   + 2*x*exp(-x)+ 3*x +3*x.*x + tan(x);

funtest3.m     5.3问题所用到的函数
function f=funtest3(x,y)
% f = -y + cos(2*x) - 2*sin(2*x) + 2*x*exp(-x);
% [x,y]=adams_wjl('funtest3',0,2,1,0.2);
% s= dsolve('Dy=-y +cos(2*x) - 2*sin(2*x) + 2*x*exp(-x)','y(0)=1','x')
% s =
% cos(2*x)+exp(-x)*x^2

%hold on
% plot(x,cos(2*x)+exp(-x).*x.^2,'*r')
f = -y + cos(2*x) - 2*sin(2*x) + 2*x*exp(-x);

%error is y-(cos(2*x)+exp(-x).*x.^2) ,not  y-cos(2*x)+exp(-x).*x.^2
% but y-(cos(2*x)-exp(-x).*x.^2)
%remember the  sign of the expression  , or the  ()

na52test.m     数值分析5.2的求解例子
%function na52test

tic
[x1,y1]=rk_wjl('funtest2',0,2,1,0.001);
t1=toc

tic
[x2,y2]=adams_wjl1('funtest2',0,2,1,0.001);
t2=toc

na53test.m     数值分析5.3的求解例子
%function na53test
[x1,y1]=rk_wjl('funtest3',0,2,1,0.1);
[x2,y2]=adams_wjl1('funtest3',0,2,1,0.1);
ye=feval('exact53',x1);
err1=y1-ye;
%err1=abs(y1-ye);
%err2=y2-ye;
err2=abs(y2-ye);
figure(1)
fprintf('%5.3f\n',x1(1:10))
fprintf('%+2.1e\n',err1(1:10))
fprintf('%+2.1e\n',err2(1:10))
plot(x1,err1,x1,err2)
legend('Runge-Kutta','Adams')
xlabel('x')
ylabel('y')
title('step=0.1,Error of problem3 with Runge-Kutta,Adams')


[x1,y1]=rk_wjl('funtest3',0,2,1,0.01);
[x2,y2]=adams_wjl1('funtest3',0,2,1,0.01);
ye=feval('exact53',x1);
fprintf('%25.15f\n',y1(1:10))
fprintf('%25.15f\n',y2(1:10))
fprintf('%25.15f\n',x1(1:10))
fprintf('%25.15f\n',ye(1:10))

err1=y1-ye;
%err1=abs(y1-ye);
err2=y2-ye;
%err2=abs(y2-ye);
figure(2)
fprintf('%+2.1e\n',err1(1:10))
fprintf('%+2.1e\n',err2(1:10))
fprintf('%5.3f\n',x1(1:10))
plot(x1,err1,x1,err2)
legend('Runge-Kutta','Adams')
xlabel('x')
ylabel('y')
title('step=0.01,Error of problem3 with Runge-Kutta,Adams')

[x1,y1]=rk_wjl('funtest3',0,2,1,0.001);
[x2,y2]=adams_wjl1('funtest3',0,2,1,0.001);
ye=feval('exact53',x1);
err1=y1-ye;
%err1=abs(y1-ye);
%err2=y2-ye;
err2=abs(y2-ye);
figure(3)
fprintf('%+2.1e\n',err1(1:10))
fprintf('%+2.1e\n',err2(1:10))
fprintf('%5.3f\n',x1(1:10))
plot(x1,err1,x1,err2)
legend('Runge-Kutta','Adams')
xlabel('x')
ylabel('y')
title('step=0.001,Error of problem3 with Runge-Kutta,Adams')



[x1,y1]=rk_wjl('funtest3',0,30,1,0.1);
[x2,y2]=adams_wjl1('funtest3',0,30,1,0.1);
ye=feval('exact53',x1);
err1=y1-ye;
%err1=abs(y1-ye);
err2=y2-ye;
%err2=abs(y2-ye);
figure(4)
fprintf('%5.3f\n',x1(1:10))
fprintf('%+2.1e\n',err1(1:10))
fprintf('%+2.1e\n',err2(1:10))
plot(x1,err1,x1,err2)

legend('Runge-Kutta','Adams')
xlabel('x')
ylabel('y')
title('step=0.1,Range=0..30,Error of problem3 with Runge-Kutta,Adams')

figure(5)
plot(x1,y1,x2,y2,x1,ye)
legend('Runge-Kutta','Adams','exact')
xlabel('x')
ylabel('y')
title('step=0.1,Range=0..30,problem3 with Runge-Kutta,Adams')

[x1,y1]=rk_wjl('funtest3',0,30,1,0.01);
[x2,y2]=adams_wjl1('funtest3',0,30,1,0.01);
ye=feval('exact53',x1);
err1=y1-ye;
%err1=abs(y1-ye);
err2=y2-ye;
%err2=abs(y2-ye);
figure(6)
fprintf('%5.3f\n',x1(1:10))
fprintf('%+2.1e\n',err1(1:10))
fprintf('%+2.1e\n',err2(1:10))
plot(x1,err1,x1,err2)

legend('Runge-Kutta','Adams')
xlabel('x')
ylabel('y')
title('step=0.01,Range=0..30,Error of problem3 with Runge-Kutta,Adams')

figure(7)
plot(x1,y1,x2,y2,x1,ye)
legend('Runge-Kutta','Adams','exact')
xlabel('x')
ylabel('y')
title('step=0.01,Range=0..30,problem3 with Runge-Kutta,Adams')

%time compare
tic;
[x1,y1]=rk_wjl('funtest3',0,2,1,0.001);
t1=toc;

tic;
[x2,y2]=adams_wjl1('funtest3',0,2,1,0.001);
t2=toc;
t1,t2

rk51test.m     数值分析5.1的求解例子
%function rk51test
%runge-kutta problem 51 test
[x1,y1]=rk_wjl('funtest1',0,1,1,0.01,50);
[x2,y2]=rk_wjl('funtest1',0,1,1,0.01,1);
[x3,y3]=rk_wjl('funtest1',0,1,1,0.01,-1);
[x4,y4]=rk_wjl('funtest1',0,1,1,0.01,-50);

y2e=exact51(x2,1);
err2=y2-y2e;
% fprintf('%5.3e\n',err2)
%plot(x2,err2)
y3e=exact51(x3,-1);
err3=y3-y3e;
%  fprintf('%5.3e\n',err3)
y4e=exact51(x4,50);
err4=y4-y4e;
%  fprintf('%5.3e\n',err4)
%  plot(x4,y4e)

figure(1)
plot(x1,y1,'--',x2,y2,'-.',x3,y3,'-..',x4,y4,'-');
legend(['\alpha=50 ';'\alpha=1  ';'\alpha=-1 ';'\alpha=-50']);
title('dy/dt=\alpha*y-\alpha*x+1,\alpha=50,1,-1,-50');
xlabel('x')
ylabel('y')

figure(2)
plot(x1,y1,'--',x2,y2,'-.',x3,y3,'-..',x4,y4,'-');
legend(['\alpha=50 ';'\alpha=1  ';'\alpha=-1 ';'\alpha=-50']);
title('dy/dt=\alpha*y-\alpha*x+1,\alpha=50,1,-1,-50');
ylim([-1,30]);
xlabel('x')
ylabel('y')

%to show the error
y1e=exact51(x1,50);
figure(3)
plot(x1,y1,x1,y1e,'*',x1,y1-y1e)
%ylim([-100,100])
title('The exact solution,\alpha=50,and the error');
legend('\alpha=50','exact solution','error');
xlabel('x')
ylabel('y')

figure(4)
plot(x1,y1,x1,y1e,'.',x1,y1-y1e)
ylim([-100,100])
legend('\alpha=50','exact solution','error');
title('The exact solution,\alpha=50,and the error,Zoom in');
fprintf('%+25.15f\n',y1-y1e)
xlabel('x')
ylabel('y')

figure(5)
error=y1-y1e;
plot(x1(1:20),error(1:20),'*')
title('Error')
xlabel('x')
ylabel('y')


figure(6)
errorbar(x1,y1e,error,':')
title('errorbar')
xlabel('x')
ylabel('y')

[x5,y5]=rk_wjl('funtest1',0,1,1,0.01,-50);
[x6,y6]=rk_wjl('funtest1',0,1,1,0.08,-50);
figure(7)
plot(x5,y5,x6,y6)
title('\alpha=-50,step with 0.01 and 0.08')
ylim([-1,50])
legend('step 0.01','step 0.08')
xlabel('x')
ylabel('y')

index=1:8:80;
%y5(index)
%y6
%fprintf('%25.15f,%25.15f\n',y5(index),y6)
fprintf('%25.15f\n',y5(index))
fprintf('%25.15f\n',y6(1:10))
%fprintf('%25.15f,%25.15f\n',y5(index),y6(1:10))
yy5=exact51(x5,-50);
fprintf('%25.15f\n',yy5(index))
fprintf('%25.15f\n',x5(index))

vanderpolequ.m VDP方程的函数
function dydt = vanderpolequ(t,y,Mu)
%VANDERPOLEQU Defines the van der Pol equation for VDPTEST.

dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];

vdptest.m      VDP求解的例子
% function vdptest
tspan = [0, 20];
y0 = [2; 0];
Mu = 1;
[t,y] = ode45(@vanderpolequ, tspan, y0,[],Mu);

% Plot of the solution
figure(1)
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')

figure(2)
plot(y(:,1),y(:,2))
xlabel('diff(y)')
ylabel('diff(y,2)')
title('phase figure')

%stiff problem
tspan = [0, 3000];
y0 = [2; 0];
Mu = 1000;
[t,y] = ode15s(@vanderpolequ, tspan, y0,[],Mu);


figure(3)
plot(t,y(:,1))
title('van der Pol Equation, \mu = 1000,solved with ode15s')
axis([0 3000 -3 3])
xlabel('t')
ylabel('solution y')

⌨️ 快捷键说明

复制代码 Ctrl + C
搜索代码 Ctrl + F
全屏模式 F11
切换主题 Ctrl + Shift + D
显示快捷键 ?
增大字号 Ctrl + =
减小字号 Ctrl + -