Como en ejercicios anteriores, para encontrar el sistema diferencial de primer orden equivalente al p.v.i. dado haremos $$x_1(t)=y(t) \hspace{.5cm}, \hspace{.5cm} x_2(t)=y'(t)$$ Escribe el sistema y pulsa en 'Ver'.

Resolución del segundo apartado

Vamos a pensar ahora cómo estimaríamos la solución de ese sistema para $r=1$ y $a=1$ utilizando 'ode45' para $t$ variando desde 0 hasta 10, con paso $h=0.1$. Recuerda que los argumentos de entrada obligatorios de 'ode45' para la aproximación de la solución de un sistema diferencial de dimensión 'm' (número de incógnitas o de ecuaciones) son

• la función $f$ definida previamente con la '@', que será un vector columna de 'm' funciones (las variables de esas funciones son: la variable independiente y las 'm' variables incógnitas);
• el vector fila de valores de $t$ donde quiera aproximarse la solución, que debe empezar en $t_0$
• el vector fila formado por los datos iniciales (valores de las incógnitas en $t_0$)

Escribe el código adecuado en este caso y pulsa en 'Ver'.

f=@(t,x) [x(2) ; -x(1)-x(1)^3];
[t,solsis]=ode45(f,0:.1:10,[1,0]);

Prepara ahora una función, 'auxDuffing', con dos variables de entrada ($r$ y $a$) y una de salida que sea la estimación de $y$ en los mismos valores de $t$ que antes. ¿Cuál de éstas es correcta?

function auxDuffing(r,a)
f=@(t,x) [x(2) ; -x(1)-r*x(1)^3];
[t,solsis]=ode45(f,0:.1:10,[a,0]);
end

function solsis=auxDuffing(r,a)
f=@(t,x) [x(2) ; -x(1)-r*x(1)^3];
[t,solsis]=ode45(f,0:.1:10,[a,0]);
end

En el fichero 'auxDuffing.m' escribimos

function sol=auxDuffing(r,a)
f=@(t,x) [x(2) ; -x(1)-r*x(1)^3];
[t,solsis]=ode45(f,0:.1:10,[a,0]);
sol=solsis(:,1);
end

y mejor si añadimos un comentario sobre lo que realiza la función:

function sol=auxDuffing(r,a)
% utiliza ode45 para estimar la solución del p.v.i y''+y+ry^3=0, y(0)=a,
% y'(0)=0; sol es el vector columna 101x1 con esa estimación
f=@(t,x) [x(2) ; -x(1)-r*x(1)^3];
[t,solsis]=ode45(f,0:.1:10,[a,0]);
sol=solsis(:,1);
end

Cuando en la ventana de comandos ejecutemos

>> sol=auxDuffing(1,1);
>> plot(0:.1:10,sol)

en la ventana de figuras se representará la estimación del problema $$y''+y+y^3=0 \ ,\ \ y(0)=1 \ ,\ \ y'(0)=0$$ para $t$ entre 0 y 10:

Continuar

Resolución del tercer apartado

En este apartado definiremos una nueva función que utilizará la función definida antes, 'auxDuffing', para dibujar las estimaciones correspondientes a varios valores de $r$ y uno fijado de $a$. Las variables de entrada de esta nueva función son por tanto el vector formado por los valores de $r$ y el valor de $a$; no tiene variable de salida y su efecto es el de dibujar tantas aproximaciones como valores de $r$ haya en el vector de entrada. Seguiremos los siguientes pasos

• definir la función con sus variables de entrada;
• definir una matriz de ceros; en ella guardaremos por columnas las estimaciones de $y(t)$;
• hacer un ciclo 'for' que para cada componente del vector de valores de $r$
  • calcule la estimación correspondiente ejecutando la función 'auxDuffing' con las variables de entrada apropiadas
  • guarde esa estimación en una columna de la matriz definida antes;
• generar el vector de valores de $t$;
• dibujar cada columna de la matriz de estimaciones respecto del vector de valores de $t$;
• cerrar la función.

Sigue esos pasos para escribir la función y pulsa en 'Ver'.

function ecuDuffing(vr,a)
% para los valores del parámetro r contenidos en vr se estima con ode45 la
% solución de y''+y+ry^3=0, y(0)=a, y'(0)=0 para los valores de t=0:.1:10
% las estimaciones se guardan como columnas de est, que es una matriz 101xlongitud(vr)
est=zeros(101,length(vr));
    for k=1:length(vr)
        sol=auxDuffing(vr(k),a);
        est(:,k)=sol;
    end
t=0:.1:10;
plot(t,est)
end

ecuDuffing(1:3,1)

Resolución del cuarto apartado

El trabajo ahora es paralelo al realizado en el apartado anterior. En este caso hemos de preparar una función que para un mismo valor de $r$, dibuje las estimaciones para distintos valores de la condición inicial $y(0)$. Habiendo resuelto el apartado anterior, éste es muy sencillo. Escribe la nueva función en un nuevo fichero y pulsa en 'Ver'.

function ecuDuffing2(r,va)
% para los valores 'a' contenidos en va se estima con ode45 la
% solución de y''+y+ry^3=0, y(0)=a, y'(0)=0 para los valores de t=0:.1:10
% las estimaciones se guardan como columnas de est, que es una matriz 101xlongitud(va)
est=zeros(101,length(va));
    for k=1:length(va)
        sol=auxDuffing(r,va(k));
        est(:,k)=sol;
    end
t=0:.1:10;
plot(t,est)
end

ecuDuffing2(1,1:4)