Solving ODEs in Matlab without ODE45 or Built-In Functions

[arizonian]

I am having trouble with some Matlab code that I am writing. I am attempting to solve ODE writing my own code, per the Professor. That means I cannot use ODE45 or other built-in functions.

My problem: Oiler, a play on Euler, calls another function, FX, that contains the two equations that I am attempting to solve. Oiler is called from the command line by the operator (me).

Command line:

[X,Y] = oiler('bloodtest',2,3,3,4)

Function # 1:

function [rhs] = oiler(FX,a,b,n,h)
r=[a;b]
s=ones(1,n)
t=ones(1,n)
rhs=[s;t]
rhs(:,1)=feval(FX,r);
for i=2:n;
rhs(:,i)=rhs(:,i-1)+feval(FX,rhs(:,i-1))*h
i=i+1;
end

Function # 2:

function [dx,dy] = bloodtest(x,y)
[dx,dy]=[x.^2;2.*y];


I get error flags saying 'too many outputs' for function #2, but I have not yet found a way to call the second function using a vector as the argument.

Any help will be appreciated

Bill

 

[Clausius2]

Function # 2:

function [dx,dy] = bloodtest(x,y)
[dx,dy]=[x.^2;2.*y];<----LOOK HERE


I get error flags saying 'too many outputs' for function #2, but I have not yet found a way to call the second function using a vector as the argument.

Hi Bill, I discovered what is reddened in the quoted text. Change the ";" in the second term by ",". Maybe that is giving you some trouble.

 

[gerben]

I think the problem is here:

rhs=[s;t]
rhs(:,1)=feval(FX,r);

in the first line you make the variable rhs a matrix with two rows and 3 columns, namely:
$$\left(\begin{array}{ccc} 1&1&1\\ 1&1&1\end{array} \right)$$
but in the next line you say: assign the result of feval(FX,r) to "all rows of the first column" of rhs. But that is a matrix with two rows and one column, while FX (bloodtest) will return a matrix with one row and two columns.

I think you should change FX (bloodtest) to:

Function # 2:

function [dx; dy] = bloodtest(x,y)
[dx; dy]=[x.^2;2.*y];

 

[arizonian]

Got past the first problem, but I am stuck on another problem. The indices don't seem to be working like I think they should.

The equations of interest from the homework assigned:
dy/dx = -2y+5*exp(-x) % (y1)
dz/dx = (-y*z^2)/2 % (y2)

I am incrementing by 0.2. They are coupled, so vector algebra comes into play. Now, on to the code that I am writing.

Command line:

EDU>> e=oiler('P25_11',0,2,4,5,.2)

Function # 1 %%%%%%%%%%%%%%%%%%%%%%%%%

function [rhs] = oiler(FX,xo,y1,y2,n,h)

rhs=ones(2,n);
rhs(:,1)=[y1;y2];
x=xo+h;
% rhs(:,1)=feval(FX,x,rhs(1,1),r(2,1));
% rhs(:,1)=feval(FX,xo,yo,zo);
for k=2:n;
rhs(:,k)=rhs(:,k-1)+(feval(FX,x,rhs(1,k-1),rhs(2,k-1)))*h;
x=x+h;
end

Function # 2 %%%%%%%%%%%%%%%%%%%%%%%

function[dy1y2]=P25_11(x,y1,y2)
dy1y2=[-2*y1+5*exp(-x);(-y1*y2^2)/2];


and the response is:

e =

2.0000 2.0187 1.8816 1.6777 1.4560
4.0000 0.8000 0.6708 0.5861 0.5285


Maybe the response is correct for what I am trying to do, but computing by hand, the first term on the bottom line should be -16. That means the entire bottom row is wrong.

Thanks again

Bill

 

[arizonian]

Maybe my understanding of differential equations with boundaries is suspect. The code does run, the incrementation is taking place like it should. I did compute the second term and it is correct. The third term is correct, and the fourth, and the fifth.

I have egg on my face.

Bill