octave-maintainers
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

## a deval function for octave

 From: philippe Subject: a deval function for octave Date: Fri, 24 May 2019 11:55:14 +0200 User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:60.0) Gecko/20100101 Thunderbird/60.6.1

```Hi to all

the syntax to solve a differential equation in matlab has recently
changed and became incompatible with the octave one ! Actually you
should do :

%% solving y''(t)+y(t)=0  y(0)=y'(0)=1 => y(t)=cos(t)+sin(t)
y0=[1;1];T=40;  %%
t=[0:0.1:T];  %%  you choose times t(k)  where solution is evaluated
sol=ode45(@(t,y) [y(2);-y(1)],[0,T],y0);   %%  solving (E)
tt=sol.x;y=sol.y;   %% retrieve tt and y(tt) vectors
clf ;  %% plotting numerical and exact solutions
plot(tt,y(1,:),'xr',t,cos(t)+sin(t),'-b')

but the problem is that you can’t choose times tt computed with ode45 ,
if you want more points you need to do an interpolation between those
computed by ode45 … in matlab this is done by deval function :

y0=[1;1];T=40;%%
t=[0:0.1:T];%%  you choose times t(k)  where solution is evaluated
sol=ode45(@(t,y) [y(2);-y(1)],[0,T],y0);
y=deval(t,sol);
clf ; %% plotting numerical and exact solutions
plot(t,y(1,:),'xr',t,cos(t)+sin(t),'-b')

but deval has not yet been implemented in octave !!! So I’ve tried to
implement my own one (see below), can it be useful to share it here?

sincerely yours ,
Philippe

function y=deval(t,sol)
%% y=deval(t,sol)
%% sol= solution of ode y'(t)=f(t,y(t))  y(0)=y0
%%       obtained with sol=ode45(@f,[t0,T],y0)
%%   t=  vector of times [t(1) ...t(n)]
%%   y= vector of solution values [y(t(1)) ... y(t(n))]
%%      interpolated  from sol data

Y=sol.y;
x=sol.x;
n=length(t);
l=size(Y,1);
y=zeros(l,n);
K=5;
for k=1:l
x_tmp=x(1:3);
ind=find((t<x(2)));
P=polyfit(x_tmp,Y(k,1:3),K);
new_y=polyval(P,t(ind)) ;
for p=2:length(x)-2
x_tmp=x(p-1:p+2);
ind=find((t>=x(p))&(t<x(p+1)));
P=polyfit(x_tmp,Y(k,p-1:p+2),K);
new_y=[new_y polyval(P,t(ind))];
end
x_tmp=x(end-2:end);
ind=find((t>=x(end-1)));
P=polyfit(x_tmp,Y(k,end-2:end),K);
new_y=[new_y polyval(P,t(ind))];
y(k,:)=new_y;
end

```

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]