solving a linear initial value problem in MATLAB

I have to solve an equation of the form: u''' = au'' + bu' + cu, u(t0) = u0, u'(t0) = ud0, u''(t0) = uddo. I have to write a function file with input u_vec: = [u0 ud0 udd0]' abc_vec: =[a b c]', t0 = initial time, tf = final time, n= number of time steps. And output: u - is a row vector of size n+1 with the first entry of u being u0 and with the (i+1) the entry of u being ui =u(ti) where i =1:100. And t = is a row vector of length n+1 containing t0....tn.

I am just stuck at where to even start with this.

Also I have to write a script file which will produce a graph of u against t. I have not a clue about how to.

I would be extremely thankful for any help.

Thank you.

Answer 1

[Jan]

At first you have to convert the system to degree 1 - this is implied in "u_vec = [u, ud, udd]" already. Then you need an integrator. Try ODE45 and check the examples include in the help and doc.

Answer 2

[Matt Tearle]

1. Convert the third-order equation into a (3-D) first-order system
2. Write a function that defines the rate equations. The system should now look like z' = f(t,z), where z = [u, u', u'']. Write your function file so that it takes t and a vector z as inputs, and returns a vector dz that represents the derivatives ([u', u'', u'''])
3. Use ode45 to solve the system, giving a function handle to you function from step 2 as input
4. Extract the appropriate component of the solution returned by ode45 to plot as a function of t (also returned by ode45).

[Lewis Watson]

thanks to both of you, so far I have this

function [u t] = test(u_vec, abc_vec, t0, tf, n)

u_vec=[u_vec]';

abc_vec=[abc_vec]';

h = tf-to/n;

x=abc_vec(1,:);

y=abc_vec(2,:);

z=abc_vec(3,:);

M = [0 1 0; 0 0 1;z y x];

for i = 1:n

t(i) = t0 + i*h;

end

t = [t0, t(i)]

am I going the right way?

Answer 3

[Matt Tearle]

Oh, sorry, it looks like Jan & I both misread/misinterpreted your questions. You're supposed to write your own ode integrator? OK, I assume, then, that this is some kind of homework problem. If not, the short answer is: don't -- use ode45 instead.

So your function above is a reasonable start. A few things, though:

1. You don't need a loop to calculate t. Use the range (:) operator or the linspace function.
2. MATLAB uses standard mathematical order of operations, so h = tf-t0/n will calculate t0/n first, then subtract from tf
3. You don't need to extract x, y, and z individually to make M. You can use fliplr and concatenation with abc_vec directly.

Next you need to write a loop to do the actual integration, using a method of your choosing. Make sure you preallocate space for the result u before you do so; then in the loop, calculate u(k) based on u(k-1). (This is much more efficient that calculating the kth u and appending it to the ones you have so far.)

[Lewis Watson]

sorry, that was my fault on the way I asked the question, thank you both for your help and time.

Answer 4

[Jia]

Hi http://www.mathworks.com/matlabcentral/answers/contributors/2764975-lewis-watson,

Can i know what is the final matlab code you get from that question?