Singular Boundary Value Problems
>function deriv(x,y) := [y[2],g3(x)-g2(x,y[1])*y[1]-g1(x,y[1])*y[2]]
>function g1(x,y) := 0
>function g2(x,y) := -phi*y^n/y
>function g3(x) := 0
>function f0(ydguess) := lsoda("deriv",[0,1],[0.1,ydguess])[1,2]
>phi=4; n=1;
>f0(0), f0(1)
0.376219569124
2.18964977308
>res=secant("f0",0,2,y=1)
0.343978185375
>f0(res)
1
>t=0:0.1:1; y=lsoda("deriv",t,[0.1,res]); ...
>plot2d(t,y[1],color=10,points=1,style="o"); ...
>insimg(25,"");

>function g1(x,y) &= 0
0
>function g2(x,y) &= -phi*y^n/y
n - 1
- phi y
>function g3(x) &= 0
0
>equ &= 'diff(y,x,2) + g1(x,y)*'diff(y,x) + g2(x,y)*y = g3(x)
2
d y n
--- - phi y = 0
2
dx
>res &= ode2(subst([n=1,phi=4],equ),y,x)
2 x - 2 x
y = %k1 E + %k2 E
>res1 &= rhs(bc2(res,x=0,y=1/10,x=1,y=1))
2 2 x 4 2 - 2 x
(10 E - 1) E (E - 10 E ) E
---------------- + -------------------
4 4
10 E - 10 10 E - 10
>x=0:0.01:1; plot2d(x,res1(x),add=true,color=13); ...
>insimg(25,"");

>function f(x,y) := [ y[2], 2*y[2]/x-4*y[1] ]
>function g(der) := lsoda("f",[1,0],[1,der])[1,-1]
>g(1)
0.236949825926
>der1 = secant("g",2.5)
2.0884291999
>x=1:-0.01:0; y=lsoda("f",x,[1,der1]); plot2d(x,y[1]); insimg;

>x=1:-0.01:-1; y=lsoda("f",x,[1,3]); plot2d(x,y[1]); ...
>x=1:-0.01:-1; y=lsoda("f",x,[1,1]); plot2d(x,y[1],add=true); ...
>insimg;

>function f(x,y) := [ y[2], 2*y[2]/x-4 ]
>function g(der) := lsoda("f",[1,0],[1,der])[1,-1]
>der1 = secant("g",2.5)
1.00000000001
>x=1:-0.01:0; y=lsoda("f",x,[1,der1]); plot2d(x,y[1]); insimg;

>sol &= ode2('diff(y,x,2)=2*'diff(y,x)/x-4,y,x)
3 2 %k1
y = %k2 x + 2 x - ---
3
>function y(x,%k1,%k2) &= rhs(sol)
3 2 %k1
%k2 x + 2 x - ---
3
>&diff(y(x,a,b),x,2) = expand(-2*diff(y(x,a,b),x)/x-4)
6 b x + 4 = - 6 b x - 12
>&solve(y(1,a,b)=1,[a,b]); ...
> function y1(x,c) &= y(x,a,b) with %[1] with %rnum_list[1]=c
3 2 - 3 c - 3
c x + 2 x + ---------
3
>&expand(y1(1,c))
1
>&diff(y1(x,c),x) with x=1
3 c + 4
>plot2d("y1(x,-4/3)",0,1); insimg;

>sol &= ode2('diff(y,x,2)=-2*'diff(y,x)/x-4,y,x)
2
2 x %k2
y = - ---- + --- + %k1
3 x
>&sqrt(x); &diff(%,x,2) + 1/2 *diff(%,x)/x
0
>function f(x,y) := [ y[2], -1/2*y[2]/x ]
>x=1:-0.01:0.01; ...
>y1=lsoda("f",x,[1,1/4])[1]; ...
>y2=lsoda("f",x,[1,3/4])[1]; ...
>plot2d(x,y1_y2); insimg;

>function h(der) := lsoda("f",[1,1e-4],[1,der])[1,2];
>der0:=bisect("h",1/4,3/4)
0.50505050494
>function f(x,y) := [ y[2], -2*y[2]/x + 4*y[1] ]
>x=1:-0.01:0.01; ...
>y1=lsoda("f",x,[1,1.2])[1]; ...
>y2=lsoda("f",x,[1,1])[1]; ...
>plot2d(x,y1_y2); insimg;

>function h(der) := lsoda("f",[1,0.01],[1,der])[1,2];
>der0:=bisect("h",1,1.2)
1.07773431527
>y=lsoda("f",x,[1,der0])[1]; ...
>plot2d(x,y); insimg;

>function f(x,a,b) &= a*exp(-2*x)/x + b*exp(2*x)/x
2 x - 2 x
b E a E
------ + --------
x x
>sol &= solve([f(1,a,b)=1,f(s,a,b)=0],[a,b])
4 s + 2 2
E E
[[a = - ---------, b = ---------]]
4 4 s 4 4 s
E - E E - E
>function fsol(x,s) &= factor(f(x,a,b) with sol[1])
2 - 2 x x s x s 2 x 2 s
E (E - E ) (E + E ) (E + E )
- ------------------------------------------
s s 2 s 2
(E - E) (E + E) (E + E ) x
>&ratsimp(diff(fsol(x,s),x,2)+2*diff(fsol(x,s),x)/x-4*fsol(x,s))
0
>plot2d("fsol(x,0.01)",color=red,add=true); insimg;

>plot2d("fsol(x,0)",0,1); insimg;

>&ratsimp(fsol(x,0))
- 2 x 4 x + 2 2
E (E - E )
----------------------
4
(E - 1) x
>function fd(x,a,b) &= diff(f(x,a,b),x)
2 x 2 x - 2 x - 2 x
2 b E b E 2 a E a E
-------- - ------ - ---------- - --------
x 2 x 2
x x
>sol &= solve([f(1,a,b)=1,fd(s,a,b)=0],[a,b])
4 s + 2
(2 s - 1) E
[[a = -----------------------------,
4 s 4
(2 s - 1) E + E (2 s + 1)
2
E (2 s + 1)
b = -----------------------------]]
4 s 4
(2 s - 1) E + E (2 s + 1)
>function fsol(x,s) &= factor(f(x,a,b) with sol[1])
2 - 2 x 4 x 4 x 4 s 4 s
E (2 s E + E + 2 s E - E )
--------------------------------------------
4 s 4 s 4 4
(2 s E - E + 2 E s + E ) x
>plot2d("fsol(x,0.1)",0,1); insimg;

>&fsol(x,0)
2 - 2 x 4 x
E (E - 1)
-------------------
4
(E - 1) x
Examples Homepage