h080504.m

(有源代码)数值分析作业,本文主要包括两个部分,第一部分是常微分方程(ODE)的三个实验题,第二部分是有关的拓展讨论,包括高阶常微分的求解和边值问题的求解(BVP).文中的算法和算例都是基于Matla

hold on
plot(x1,exp(x) +x,'*r')
plot(x1,exp(x1) +x1,'*r')
cc
[x50,y50]=rk_wjl('funtest1',0,1,1,0.01,50);
yy50=exact51(x,50);
yy50=exact51(x50,50);
plot(x50,y50,x50,yy50)
[x50,y50]=rk_wjl('funtest1',0,2,1,0.01,50);
yy50=exact51(x50,50);
plot(x50,y50,x50,yy50)
plot(x50,yy50-y50)
ylim([-1,1000])
help xlswrite
cc
[x2fh001,y2fh001]=rk_wjl('funtest1',0,1,1,0.01,-2);
%-- 08-5-4 上午9:40 --%
cdtemp
cd_temp
dir
open dampmass
open odesol
odesol_w
open odesol_w
doc ode45
options = odeset('RelTol',1e-4,'AbsTol',[1e-3 2e-3 3e-3]);
odeget(options,'RelTol')
o = odeset('MaxStep',1e-6)
[x,y]=ode45('dampmass',[0 1],0.001,o);
o = odeset('MaxStep',1e-4) ;
[x,y]=ode45('dampmass',[0 1],0.001,o);
plot(x,y)
ode45(@ dampmass,[0 1],0.001);
ode45(@dampmass,[0 1],0.001);
[x,y]=ode45('dampmass',[0 1],0.001,o);
da=feval('dampmass',x,y);
[x,y]=ode45('dampmass',[0 1],0.001,o);
cc
[x,y]=ode45('dampmass',[0 1],0.001,o);
[x,y]=ode45('dampmass',[0 1],0.001);
for i=length(x)
da(i)=feval('dampmass',x(i),y(i));
end
plot(x,da)
ylim([-0.3,0.3])
y
x
y
[x,y]=ode45('dampmass',[0 1],0.001);
plot(x,y)
dydt=feval('dampmass1',y);
plot(x,dydt)
dx=dydt.*sin(75*x) + y.*cos(75*x)*75;
plot(x,dx)
plot(x,dydt)
plot(x,abs(dydt))
plot(x,log(abs(dydt)))
plot(x,log10(abs(dydt)))
[x,y]=ode45('dampmass',[0 1],0.001,o);
dir
open vdpdemo
vdpdemo
cc
help shg
open vdp1000
open vdp1
open vdpm
open vdpode
doc vdpode
cc
o = odeset('MaxStep',1e-5) ;
[x,y]=ode45('dampmass',[0 1],0.001);
size(x)
size(y)
doc  odeset
odeset(o)
o = odeset('MaxStep',1e-5,'RelTol',1e-4) ;
[x,y]=ode45('dampmass',[0 1],0.001);
o = odeset('MaxStep',1e-5) ;
[x,y]=ode45('dampmass',[0 1],0.001,o);
size(x)
plot(x,y)
plot(x,y,'*')
plot(x,y)
plot(x,y,'*')
o = odeset('MaxStep',1e-4) ;
[x,y]=ode45('dampmass',[0 1],0.001,o);
plot(x,y)
size(x)
siny=y.*sin(75*x);
plot(x,siny)
plot(x,y)
plot(x,abs(y))
plot(x,diff(real(y)))
plot(x(1:end-1),diff(real(y)))
plot(x(1:end-1),abs(diff(real(y))))
diffy=diff(real(y));
sindy=diffy.*(sin(75*x(1:end-1));
sindy=diffy.*(sin(75*x(1:end-1)));
plot(x(1:end-1),sindy)
a=[1 2 3 4];
b=[2 3 4 5];
a.*b
dydt=feval('dampmass1',y);
plot(x,y)
dydt=feval('dampmass1',y);
plot(x,y)
plot(x,dydt)
sindy=dydt.*(sin(75*x(1:end-1)));
sindy=dydt.*(sin(75*x));
plot(x,sindy)
o = odeset('MaxStep',1e-4) ;
[x,y]=ode45('dampmass',[0 1],0.001,o);
dydt=feval('dampmass1',y);
sindy=dydt.*(sin(75*x));
plot(x,sindy)
sindy=dydt.*(sin(226*x));
plot(x,sindy)
sindy=dydt.*(sin(2*pi*75*x));
plot(x,sindy)
o = odeset('MaxStep',1e-4) ;
[x,y]=ode45('dampmass',[0 1],0.001,o);
dydt=feval('dampmass1',y);
sindy=dydt.*(sin(2*pi*75*x));
plot(x,sindy)
plot(x,dydt)
plot(x,sindy)
o = odeset('MaxStep',1e-4) ;
[x,y]=ode45('dampmass',[0 1],0.001,o);
dydt=feval('dampmass1',y);
sindy=dydt.*(sin(2*pi*75*x));
cc
open vdpdemo
open odedemo
type vanderpoldemo
tspan = [0, 20];
y0 = [2; 0];
Mu = 1;
[t,y] = ode45(@vanderpoldemo, tspan, y0,[],Mu);
% Plot of the solution
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')
tspan = [0, 3000];
y0 = [2; 0];
Mu = 1000;
[t,y] = ode15s(@vanderpoldemo, tspan, y0,[],Mu);
plot(t,y(:,1))
title('van der Pol Equation, \mu = 1000')
axis([0 3000 -3 3])
xlabel('t')
ylabel('solution y')
type twoode
solinit = bvpinit([0 1 2 3 4],[1; 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
doc odedemo
doc bvp4c
solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
hold off
type twobc
open vanderpoldemo
tspan = [0, 20];
y0 = [2; 0];
Mu = 1;
[t,y] = ode45(@vanderpoldemo, tspan, y0,[],Mu);
% Plot of the solution
plot(t,y(:,1))
xlabel('t')
ylabel('solution y')
title('van der Pol Equation, \mu = 1')
figure(2)
plot(y(:,1),y(:,2))
figure(2)
plot(y(:,1),y(:,2))
xlabel('diff(y)')
ylabel('diff(y,2)')
title('phase figure')
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')
cc
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')
figure(3)
plot(t,y(:,1))
title('van der Pol Equation, \mu = 1000,solved with ode15s')
tspan = [0, 3000];
y0 = [2; 0];
Mu = 1000;
[t,y] = ode15s(@vanderpoldemo, 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')
cc
[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);
yy5=exact51(x5,-50);
fprintf('%25.15f\n',yy5(index))
fprintf('%25.15f\n',x5(index))
[x5,y5]=rk_wjl('funtest1',0,1,1,0.01,-50);
[x6,y6]=rk_wjl('funtest1',0,1,1,0.08,-50);
yy5=exact51(x5,-50);
fprintf('%25.15f\n',yy5(index))
fprintf('%25.15f\n',x5(index))
index=1:8:80;
yy5=exact51(x5,-50);
fprintf('%25.15f\n',yy5(index))
fprintf('%25.15f\n',x5(index))
fprintf('%25.15f\n',y5(index))
fprintf('%25.15f\n',y6(1:10))
y2e=exact51(x2,1);
err2=y2-y2e;
fprintf('%25.15f\n',err2)
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)
fprintf('%5.3e\n',err3)
y4e=exact51(x4,-50);
err4=y4-y4e;
fprintf('%5.3e\n',err4)
cc
[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);
y4e=exact51(x4,-50);
err4=y4-y4e;
fprintf('%5.3e\n',err4)
-50*0.01
0.01
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.
solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
hold off
solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
hold off
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
solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
hold off
cc
legend('sol','solinit')
title('BVP')
solinit = bvpinit([0 1 2 3 4],[1; 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')
title('BVP')
solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
hold off
legend('sol','solinit')
title('BVP')
solinit = bvpinit([0 1 2 3 4],[-2; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
hold off
legend('sol','solinit')
title('BVP')
figure
hold on
solinit = bvpinit([0 1 2 3 4],[-2; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
legend('sol','solinit')
title('BVP')
solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
legend('sol','solinit')
title('BVP')
solinit = bvpinit([0 1 2 3 4],[-10; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
close all
cc
hold on
solinit = bvpinit([0 1 2 3 4],[-10; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
hold on
solinit = bvpinit([0 1 2 3 4],[1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
solinit = bvpinit([0 1 2 3 4],[-1; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
solinit = bvpinit([0 1 2 3 4],[-6; 0]);
sol = bvp4c(@twoode,@twobc,solinit);
xint = linspace(0,4,50);
yint = deval(sol,xint);
plot(xint,yint(1,:),'r');
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');
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');
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');
clf
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');
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,:),'g');
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');
1e4
1e-4
3*(1e-4)
3*1e-4
3*1e-4*4
cc
notebook setup
help notebook
notebook -setup
3
notebook -setup
%-- 08-5-4 下午8:10 --%
notebook -setup
3
open adams_wjl1
open bvptest
open exact51
open exact53
open funtest1
open funtest2
open funtest3
open na52test
open na53test
open na51test
open rk51test
open vanderpolequ
cc