Solving MATLAB ODE with Matrix Coefficients: Help Requested

[christy]

Hi,

I have tried everything to figure this out in MATLAB and I can't so I was wondering if anyone could help me please.

I have a 2nd order linear system that I reduced to 1st order such that
q = dz/dt
and A.(dq/dt) = Gq(t) + f(t)

where A and G are matrices...I can't figure out how to solve this in MATLAB with having matrix coefficients. Can someone please give me some hints?

A and G are generated matrices and are constant for this problem.

Thanks

 

[LeBrad]

You can solve systems of equations and systems of systems of equations the same way.

If q is size N x 1, and z is size N x 1, then define
y = [q; z] which is size 2N x 1
Using a simple forward Euler scheme,
y[n+1] = y[n] + dt*[f1(q,z,t(n)) ; f2(q,z,t(n))]
where f2 = q[n] and f1 = inv(A)*(G*q[n] + f(tn))
G is size N x N, q is N x 1, and f is N x 1, so [f1; f2] is size 2N x 1, just like y. So all the sizes work fine. You're just solving one really big system. Then you can make some operator C which is size 1 x N to pick out the ouput you want, like z_1(t), or whatever.

Alternatively, you could keep the two equations separate and solve
q[n+1] = q[n] + dt*inv(A)*(G*q[n] + f(tn))
and
z[n+1] = z[n] + dt*(q[n])

Now you're solving two systems, each size N x 1.

If you still don't like that, you can break up q and z into N equations each.
Since q = [q_1; q_2; q_3;...;q_N], you can write it as
q_1[n+1] = q_1[n]+dt*((inv(A)*G)(1,:)*q[n]+inv(A)f_1(tn))
q_2[n+1] = q_2[n]+dt*((inv(A)*G)(2,:)*q[n]+inv(A)f_2(tn))
.
.
.
and do the same for z to get 2N individual equations to solve.

So the bottom line is, you can either stick all of your equations, each of which may be working on a vector, into one giant vector and solve that, or you can break it up into any number of equations and solve them all individually stepping through time.

 

[nkamaraju]

Dear all,
I have a equation that goes like this...
d^2/dt^2[x] + f(t) x - gamma * d/dt [x] where f(t) is a function of t. I have a numerical form of f(t), but no analytical form, so i want to pass this whole vector into my function...but I am able to pass only a single value...My program is as follows,

Main.m
clear all
close all;
clc;
s1 = input('Input File name without any *.txt extention :','s');
input_file_name = [s1 '.txt'];
Data=load(input_file_name); % column vector

time = Data(:,1)'; % Convert to row vecor
init = [0,0.5*2*pi]; % Initial condition
f =1.42E27*Data(:,2); % f(t) -- only numerical form is available

[t,y]=ode45(@eqn,time,init,[],f); % If u remove f from here and write this as %exp(-t/2E-12) inside the fun eqn.m, then, it works
plot(t,y(:,1))

eqn.m
function rk = eqn(t,y,f)
gamma= 2*pi*0.1E12;
rk = [y(2); - (f y(1)) -gamma*y(2)];

Could anyone help me out to solve this problem...

with best regards
kamaraju